DUNE PDELab (2.8)

lagrangeprism.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_LAGRANGEPRISM_HH
4#define DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGEPRISM_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 LagrangePrismLocalBasis
33 {
34 static constexpr std::size_t dim = 3;
35 public:
36 using Traits = LocalBasisTraits<D,dim,FieldVector<D,dim>,R,1,FieldVector<R,1>,FieldMatrix<R,1,dim> >;
37
40 static constexpr unsigned int size ()
41 {
42 return binomial(k+2,2u) * (k+1);
43 }
44
46 void evaluateFunction(const typename Traits::DomainType& in,
47 std::vector<typename Traits::RangeType>& out) const
48 {
49 out.resize(size());
50
51 // Specialization for zero-order case
52 if (k==0)
53 {
54 out[0] = 1;
55 return;
56 }
57
58 if (k==1)
59 {
60 out[0] = (1.0-in[0]-in[1])*(1.0-in[2]);
61 out[1] = in[0]*(1-in[2]);
62 out[2] = in[1]*(1-in[2]);
63 out[3] = in[2]*(1.0-in[0]-in[1]);
64 out[4] = in[0]*in[2];
65 out[5] = in[1]*in[2];
66
67 return;
68 }
69
70 if (k==2)
71 {
72 FieldVector<R,k+1> segmentShapeFunction;
73 segmentShapeFunction[0] = 1 + in[2] * (-3 + 2*in[2]);
74 segmentShapeFunction[1] = in[2] * (4 - 4*in[2]);
75 segmentShapeFunction[2] = in[2] * (-1 + 2*in[2]);
76
77 FieldVector<R, 6> triangleShapeFunction;
78 triangleShapeFunction[0] = 2 * (1 - in[0] - in[1]) * (0.5 - in[0] - in[1]);
79 triangleShapeFunction[1] = 2 * in[0] * (-0.5 + in[0]);
80 triangleShapeFunction[2] = 2 * in[1] * (-0.5 + in[1]);
81 triangleShapeFunction[3] = 4*in[0] * (1 - in[0] - in[1]);
82 triangleShapeFunction[4] = 4*in[1] * (1 - in[0] - in[1]);
83 triangleShapeFunction[5] = 4*in[0]*in[1];
84
85 // lower triangle:
86 out[0] = triangleShapeFunction[0] * segmentShapeFunction[0];
87 out[1] = triangleShapeFunction[1] * segmentShapeFunction[0];
88 out[2] = triangleShapeFunction[2] * segmentShapeFunction[0];
89
90 //upper triangle
91 out[3] = triangleShapeFunction[0] * segmentShapeFunction[2];
92 out[4] = triangleShapeFunction[1] * segmentShapeFunction[2];
93 out[5] = triangleShapeFunction[2] * segmentShapeFunction[2];
94
95 // vertical edges
96 out[6] = triangleShapeFunction[0] * segmentShapeFunction[1];
97 out[7] = triangleShapeFunction[1] * segmentShapeFunction[1];
98 out[8] = triangleShapeFunction[2] * segmentShapeFunction[1];
99
100 // lower triangle edges
101 out[9] = triangleShapeFunction[3] * segmentShapeFunction[0];
102 out[10] = triangleShapeFunction[4] * segmentShapeFunction[0];
103 out[11] = triangleShapeFunction[5] * segmentShapeFunction[0];
104
105 // upper triangle edges
106 out[12] = triangleShapeFunction[3] * segmentShapeFunction[2];
107 out[13] = triangleShapeFunction[4] * segmentShapeFunction[2];
108 out[14] = triangleShapeFunction[5] * segmentShapeFunction[2];
109
110 // quadrilateral sides
111 out[15] = triangleShapeFunction[3] * segmentShapeFunction[1];
112 out[16] = triangleShapeFunction[4] * segmentShapeFunction[1];
113 out[17] = triangleShapeFunction[5] * segmentShapeFunction[1];
114
115 return;
116 }
117
118 DUNE_THROW(NotImplemented, "LagrangePrismLocalBasis::evaluateFunction for order " << k);
119 }
120
126 void evaluateJacobian(const typename Traits::DomainType& in,
127 std::vector<typename Traits::JacobianType>& out) const
128 {
129 out.resize(size());
130
131 // Specialization for k==0
132 if (k==0)
133 {
134 std::fill(out[0][0].begin(), out[0][0].end(), 0);
135 return;
136 }
137
138 if (k==1)
139 {
140 out[0][0] = {in[2]-1, in[2]-1, in[0]+in[1]-1};
141 out[1][0] = {1-in[2], 0, -in[0]};
142 out[2][0] = { 0, 1-in[2], -in[1]};
143 out[3][0] = { -in[2], -in[2], 1-in[0]-in[1]};
144 out[4][0] = { in[2], 0, in[0]};
145 out[5][0] = { 0, in[2], in[1]};
146
147 return;
148 }
149
150 if (k==2)
151 {
152 // Second-order shape functions on a triangle, and the first derivatives
153 FieldVector<R, 6> triangleShapeFunction;
154 triangleShapeFunction[0] = 2 * (1 - in[0] - in[1]) * (0.5 - in[0] - in[1]);
155 triangleShapeFunction[1] = 2 * in[0] * (-0.5 + in[0]);
156 triangleShapeFunction[2] = 2 * in[1] * (-0.5 + in[1]);
157 triangleShapeFunction[3] = 4*in[0] * (1 - in[0] - in[1]);
158 triangleShapeFunction[4] = 4*in[1] * (1 - in[0] - in[1]);
159 triangleShapeFunction[5] = 4*in[0]*in[1];
160
161 std::array<std::array<R,2>,6> triangleShapeFunctionDer;
162 triangleShapeFunctionDer[0] = {-3 + 4*(in[0] + in[1]), -3 + 4*(in[0] + in[1])};
163 triangleShapeFunctionDer[1] = { -1 + 4*in[0], 0};
164 triangleShapeFunctionDer[2] = { 0, -1 + 4*in[1]};
165 triangleShapeFunctionDer[3] = { 4 - 8*in[0] - 4*in[1], -4*in[0]};
166 triangleShapeFunctionDer[4] = { -4*in[1], 4 - 4*in[0] - 8*in[1]};
167 triangleShapeFunctionDer[5] = { 4*in[1], 4*in[0]};
168
169 // Second-order shape functions on a line, and the first derivatives
170 FieldVector<R,k+1> segmentShapeFunction;
171 segmentShapeFunction[0] = 1 + in[2] * (-3 + 2*in[2]);
172 segmentShapeFunction[1] = in[2] * ( 4 - 4*in[2]);
173 segmentShapeFunction[2] = in[2] * (-1 + 2*in[2]);
174
175 FieldVector<R,k+1> segmentShapeFunctionDer;
176 segmentShapeFunctionDer[0] = -3 + 4*in[2];
177 segmentShapeFunctionDer[1] = 4 - 8*in[2];
178 segmentShapeFunctionDer[2] = -1 + 4*in[2];
179
180 // lower triangle:
181 out[0][0][0] = triangleShapeFunctionDer[0][0] * segmentShapeFunction[0];
182 out[0][0][1] = triangleShapeFunctionDer[0][1] * segmentShapeFunction[0];
183 out[0][0][2] = triangleShapeFunction[0] * segmentShapeFunctionDer[0];
184
185 out[1][0][0] = triangleShapeFunctionDer[1][0] * segmentShapeFunction[0];
186 out[1][0][1] = triangleShapeFunctionDer[1][1] * segmentShapeFunction[0];
187 out[1][0][2] = triangleShapeFunction[1] * segmentShapeFunctionDer[0];
188
189 out[2][0][0] = triangleShapeFunctionDer[2][0] * segmentShapeFunction[0];
190 out[2][0][1] = triangleShapeFunctionDer[2][1] * segmentShapeFunction[0];
191 out[2][0][2] = triangleShapeFunction[2] * segmentShapeFunctionDer[0];
192
193 //upper triangle
194 out[3][0][0] = triangleShapeFunctionDer[0][0] * segmentShapeFunction[2];
195 out[3][0][1] = triangleShapeFunctionDer[0][1] * segmentShapeFunction[2];
196 out[3][0][2] = triangleShapeFunction[0] * segmentShapeFunctionDer[2];
197
198 out[4][0][0] = triangleShapeFunctionDer[1][0] * segmentShapeFunction[2];
199 out[4][0][1] = triangleShapeFunctionDer[1][1] * segmentShapeFunction[2];
200 out[4][0][2] = triangleShapeFunction[1] * segmentShapeFunctionDer[2];
201
202 out[5][0][0] = triangleShapeFunctionDer[2][0] * segmentShapeFunction[2];
203 out[5][0][1] = triangleShapeFunctionDer[2][1] * segmentShapeFunction[2];
204 out[5][0][2] = triangleShapeFunction[2] * segmentShapeFunctionDer[2];
205
206 // vertical edges
207 out[6][0][0] = triangleShapeFunctionDer[0][0] * segmentShapeFunction[1];
208 out[6][0][1] = triangleShapeFunctionDer[0][1] * segmentShapeFunction[1];
209 out[6][0][2] = triangleShapeFunction[0] * segmentShapeFunctionDer[1];
210
211 out[7][0][0] = triangleShapeFunctionDer[1][0] * segmentShapeFunction[1];
212 out[7][0][1] = triangleShapeFunctionDer[1][1] * segmentShapeFunction[1];
213 out[7][0][2] = triangleShapeFunction[1] * segmentShapeFunctionDer[1];
214
215 out[8][0][0] = triangleShapeFunctionDer[2][0] * segmentShapeFunction[1];
216 out[8][0][1] = triangleShapeFunctionDer[2][1] * segmentShapeFunction[1];
217 out[8][0][2] = triangleShapeFunction[2] * segmentShapeFunctionDer[1];
218
219 // lower triangle edges
220 out[9][0][0] = triangleShapeFunctionDer[3][0] * segmentShapeFunction[0];
221 out[9][0][1] = triangleShapeFunctionDer[3][1] * segmentShapeFunction[0];
222 out[9][0][2] = triangleShapeFunction[3] * segmentShapeFunctionDer[0];
223
224 out[10][0][0] = triangleShapeFunctionDer[4][0] * segmentShapeFunction[0];
225 out[10][0][1] = triangleShapeFunctionDer[4][1] * segmentShapeFunction[0];
226 out[10][0][2] = triangleShapeFunction[4] * segmentShapeFunctionDer[0];
227
228 out[11][0][0] = triangleShapeFunctionDer[5][0] * segmentShapeFunction[0];
229 out[11][0][1] = triangleShapeFunctionDer[5][1] * segmentShapeFunction[0];
230 out[11][0][2] = triangleShapeFunction[5] * segmentShapeFunctionDer[0];
231
232 // upper triangle edges
233 out[12][0][0] = triangleShapeFunctionDer[3][0] * segmentShapeFunction[2];
234 out[12][0][1] = triangleShapeFunctionDer[3][1] * segmentShapeFunction[2];
235 out[12][0][2] = triangleShapeFunction[3] * segmentShapeFunctionDer[2];
236
237 out[13][0][0] = triangleShapeFunctionDer[4][0] * segmentShapeFunction[2];
238 out[13][0][1] = triangleShapeFunctionDer[4][1] * segmentShapeFunction[2];
239 out[13][0][2] = triangleShapeFunction[4] * segmentShapeFunctionDer[2];
240
241 out[14][0][0] = triangleShapeFunctionDer[5][0] * segmentShapeFunction[2];
242 out[14][0][1] = triangleShapeFunctionDer[5][1] * segmentShapeFunction[2];
243 out[14][0][2] = triangleShapeFunction[5] * segmentShapeFunctionDer[2];
244
245 // quadrilateral sides
246 out[15][0][0] = triangleShapeFunctionDer[3][0] * segmentShapeFunction[1];
247 out[15][0][1] = triangleShapeFunctionDer[3][1] * segmentShapeFunction[1];
248 out[15][0][2] = triangleShapeFunction[3] * segmentShapeFunctionDer[1];
249
250 out[16][0][0] = triangleShapeFunctionDer[4][0] * segmentShapeFunction[1];
251 out[16][0][1] = triangleShapeFunctionDer[4][1] * segmentShapeFunction[1];
252 out[16][0][2] = triangleShapeFunction[4] * segmentShapeFunctionDer[1];
253
254 out[17][0][0] = triangleShapeFunctionDer[5][0] * segmentShapeFunction[1];
255 out[17][0][1] = triangleShapeFunctionDer[5][1] * segmentShapeFunction[1];
256 out[17][0][2] = triangleShapeFunction[5] * segmentShapeFunctionDer[1];
257
258 return;
259 }
260
261 DUNE_THROW(NotImplemented, "LagrangePrismLocalBasis::evaluateJacobian for order " << k);
262 }
263
270 void partial(const std::array<unsigned int,dim>& order,
271 const typename Traits::DomainType& in,
272 std::vector<typename Traits::RangeType>& out) const
273 {
274 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
275
276 out.resize(size());
277
278 if (totalOrder == 0)
279 {
280 evaluateFunction(in, out);
281 return;
282 }
283
284 // Specialization for zero-order finite elements
285 if (k==0)
286 {
287 out[0] = 0;
288 return;
289 }
290
291 // Specialization for first-order finite elements
292 if (k==1)
293 {
294 if (totalOrder == 1)
295 {
296 auto direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
297
298 switch (direction) {
299 case 0:
300 out[0] = in[2]-1;
301 out[1] = 1-in[2];
302 out[2] = 0;
303 out[3] = -in[2];
304 out[4] = in[2];
305 out[5] = 0;
306 break;
307 case 1:
308 out[0] = in[2]-1;
309 out[1] = 0;
310 out[2] = 1-in[2];
311 out[3] = -in[2];
312 out[4] = 0;
313 out[5] = in[2];
314 break;
315 case 2:
316 out[0] = in[0]+in[1]-1;
317 out[1] = -in[0];
318 out[2] = -in[1];
319 out[3] = 1-in[0]-in[1];
320 out[4] = in[0];
321 out[5] = in[1];
322 break;
323 default:
324 DUNE_THROW(RangeError, "Component out of range.");
325 }
326 } else if (totalOrder == 2) {
327 out.resize(size());
328 if (order[0] == 1 && order[2] == 1) {
329 out[0] = 1;
330 out[1] =-1;
331 out[2] = 0;
332 out[3] =-1;
333 out[4] = 1;
334 out[5] = 0;
335 } else if (order[1] == 1 && order[2] == 1) {
336 out[0] = 1;
337 out[1] = 0;
338 out[2] =-1;
339 out[3] =-1;
340 out[4] = 0;
341 out[5] = 1;
342 } else {
343 for (std::size_t i = 0; i < size(); ++i)
344 out[i] = 0;
345 }
346 } else {
347 out.resize(size());
348 std::fill(out.begin(), out.end(), 0.0);
349 }
350
351 return;
352 }
353
354 // Specialization for second-order finite elements
355 if (k==2)
356 {
357 if (totalOrder == 1)
358 {
359 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
360 switch (direction)
361 {
362 case 0:
363 {
364 FieldVector<R,6> triangleShapeFunctionDerX;
365 triangleShapeFunctionDerX[0] = -3 + 4*(in[0] + in[1]);
366 triangleShapeFunctionDerX[1] = -1 + 4* in[0];
367 triangleShapeFunctionDerX[2] = 0;
368 triangleShapeFunctionDerX[3] = 4 - 8* in[0] - 4*in[1];
369 triangleShapeFunctionDerX[4] = -4*in[1];
370 triangleShapeFunctionDerX[5] = 4*in[1];
371
372 FieldVector<R,k+1> segmentShapeFunction;
373 segmentShapeFunction[0] = 1 + in[2] * (-3 + 2*in[2]);
374 segmentShapeFunction[1] = in[2] * ( 4 - 4*in[2]);
375 segmentShapeFunction[2] = in[2] * (-1 + 2*in[2]);
376
377 out[0] = triangleShapeFunctionDerX[0] * segmentShapeFunction[0];
378 out[1] = triangleShapeFunctionDerX[1] * segmentShapeFunction[0];
379 out[2] = triangleShapeFunctionDerX[2] * segmentShapeFunction[0];
380 out[3] = triangleShapeFunctionDerX[0] * segmentShapeFunction[2];
381 out[4] = triangleShapeFunctionDerX[1] * segmentShapeFunction[2];
382 out[5] = triangleShapeFunctionDerX[2] * segmentShapeFunction[2];
383 out[6] = triangleShapeFunctionDerX[0] * segmentShapeFunction[1];
384 out[7] = triangleShapeFunctionDerX[1] * segmentShapeFunction[1];
385 out[8] = triangleShapeFunctionDerX[2] * segmentShapeFunction[1];
386 out[9] = triangleShapeFunctionDerX[3] * segmentShapeFunction[0];
387 out[10] = triangleShapeFunctionDerX[4] * segmentShapeFunction[0];
388 out[11] = triangleShapeFunctionDerX[5] * segmentShapeFunction[0];
389 out[12] = triangleShapeFunctionDerX[3] * segmentShapeFunction[2];
390 out[13] = triangleShapeFunctionDerX[4] * segmentShapeFunction[2];
391 out[14] = triangleShapeFunctionDerX[5] * segmentShapeFunction[2];
392 out[15] = triangleShapeFunctionDerX[3] * segmentShapeFunction[1];
393 out[16] = triangleShapeFunctionDerX[4] * segmentShapeFunction[1];
394 out[17] = triangleShapeFunctionDerX[5] * segmentShapeFunction[1];
395 break;
396 }
397 case 1:
398 {
399 FieldVector<R,6> triangleShapeFunctionDerY;
400 triangleShapeFunctionDerY[0] = -3 + 4*(in[0] + in[1]);
401 triangleShapeFunctionDerY[1] = 0;
402 triangleShapeFunctionDerY[2] = -1 + 4* in[1];
403 triangleShapeFunctionDerY[3] = -4* in[0];
404 triangleShapeFunctionDerY[4] = 4 - 4* in[0] - 8*in[1];
405 triangleShapeFunctionDerY[5] = 4* in[0];
406
407 FieldVector<R,k+1> segmentShapeFunction;
408 segmentShapeFunction[0] = 1 + in[2] * (-3 + 2*in[2]);
409 segmentShapeFunction[1] = in[2] * ( 4 - 4*in[2]);
410 segmentShapeFunction[2] = in[2] * (-1 + 2*in[2]);
411
412 out[0] = triangleShapeFunctionDerY[0] * segmentShapeFunction[0];
413 out[1] = triangleShapeFunctionDerY[1] * segmentShapeFunction[0];
414 out[2] = triangleShapeFunctionDerY[2] * segmentShapeFunction[0];
415 out[3] = triangleShapeFunctionDerY[0] * segmentShapeFunction[2];
416 out[4] = triangleShapeFunctionDerY[1] * segmentShapeFunction[2];
417 out[5] = triangleShapeFunctionDerY[2] * segmentShapeFunction[2];
418 out[6] = triangleShapeFunctionDerY[0] * segmentShapeFunction[1];
419 out[7] = triangleShapeFunctionDerY[1] * segmentShapeFunction[1];
420 out[8] = triangleShapeFunctionDerY[2] * segmentShapeFunction[1];
421 out[9] = triangleShapeFunctionDerY[3] * segmentShapeFunction[0];
422 out[10] = triangleShapeFunctionDerY[4] * segmentShapeFunction[0];
423 out[11] = triangleShapeFunctionDerY[5] * segmentShapeFunction[0];
424 out[12] = triangleShapeFunctionDerY[3] * segmentShapeFunction[2];
425 out[13] = triangleShapeFunctionDerY[4] * segmentShapeFunction[2];
426 out[14] = triangleShapeFunctionDerY[5] * segmentShapeFunction[2];
427 out[15] = triangleShapeFunctionDerY[3] * segmentShapeFunction[1];
428 out[16] = triangleShapeFunctionDerY[4] * segmentShapeFunction[1];
429 out[17] = triangleShapeFunctionDerY[5] * segmentShapeFunction[1];
430 break;
431 }
432 case 2:
433 {
434 FieldVector<R, 6> triangleShapeFunction;
435 triangleShapeFunction[0] = 2 * (1 - in[0] - in[1]) * (0.5 - in[0] - in[1]);
436 triangleShapeFunction[1] = 2 * in[0] * (-0.5 + in[0]);
437 triangleShapeFunction[2] = 2 * in[1] * (-0.5 + in[1]);
438 triangleShapeFunction[3] = 4*in[0] * (1 - in[0] - in[1]);
439 triangleShapeFunction[4] = 4*in[1] * (1 - in[0] - in[1]);
440 triangleShapeFunction[5] = 4*in[0]*in[1];
441
442 FieldVector<R,k+1> segmentShapeFunctionDer;
443 segmentShapeFunctionDer[0] = -3 + 4*in[2];
444 segmentShapeFunctionDer[1] = 4 - 8*in[2];
445 segmentShapeFunctionDer[2] = -1 + 4*in[2];
446
447 out[0] = triangleShapeFunction[0] * segmentShapeFunctionDer[0];
448 out[1] = triangleShapeFunction[1] * segmentShapeFunctionDer[0];
449 out[2] = triangleShapeFunction[2] * segmentShapeFunctionDer[0];
450 out[3] = triangleShapeFunction[0] * segmentShapeFunctionDer[2];
451 out[4] = triangleShapeFunction[1] * segmentShapeFunctionDer[2];
452 out[5] = triangleShapeFunction[2] * segmentShapeFunctionDer[2];
453 out[6] = triangleShapeFunction[0] * segmentShapeFunctionDer[1];
454 out[7] = triangleShapeFunction[1] * segmentShapeFunctionDer[1];
455 out[8] = triangleShapeFunction[2] * segmentShapeFunctionDer[1];
456 out[9] = triangleShapeFunction[3] * segmentShapeFunctionDer[0];
457 out[10] = triangleShapeFunction[4] * segmentShapeFunctionDer[0];
458 out[11] = triangleShapeFunction[5] * segmentShapeFunctionDer[0];
459 out[12] = triangleShapeFunction[3] * segmentShapeFunctionDer[2];
460 out[13] = triangleShapeFunction[4] * segmentShapeFunctionDer[2];
461 out[14] = triangleShapeFunction[5] * segmentShapeFunctionDer[2];
462 out[15] = triangleShapeFunction[3] * segmentShapeFunctionDer[1];
463 out[16] = triangleShapeFunction[4] * segmentShapeFunctionDer[1];
464 out[17] = triangleShapeFunction[5] * segmentShapeFunctionDer[1];
465 break;
466 }
467 default:
468 DUNE_THROW(RangeError, "Component out of range.");
469 }
470 } else {
471 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
472 }
473
474 return;
475 }
476
477 DUNE_THROW(NotImplemented, "LagrangePrismLocalBasis::partial not implemented for order " << k);
478 }
479
481 static constexpr unsigned int order ()
482 {
483 return k;
484 }
485 };
486
491 template<unsigned int k>
492 class LagrangePrismLocalCoefficients
493 {
494 public:
496 LagrangePrismLocalCoefficients ()
497 : localKeys_(size())
498 {
499 if (k==0)
500 {
501 localKeys_[0] = LocalKey(0,0,0);
502 return;
503 }
504
505 if (k==1)
506 {
507 for (std::size_t i=0; i<size(); i++)
508 localKeys_[i] = LocalKey(i,3,0);
509 return;
510 }
511
512 if (k==2)
513 {
514 // Vertex shape functions
515 localKeys_[0] = LocalKey(0,3,0);
516 localKeys_[1] = LocalKey(1,3,0);
517 localKeys_[2] = LocalKey(2,3,0);
518 localKeys_[3] = LocalKey(3,3,0);
519 localKeys_[4] = LocalKey(4,3,0);
520 localKeys_[5] = LocalKey(5,3,0);
521
522 // Edge shape functions
523 localKeys_[6] = LocalKey(0,2,0);
524 localKeys_[7] = LocalKey(1,2,0);
525 localKeys_[8] = LocalKey(2,2,0);
526 localKeys_[9] = LocalKey(3,2,0);
527 localKeys_[10] = LocalKey(4,2,0);
528 localKeys_[11] = LocalKey(5,2,0);
529 localKeys_[12] = LocalKey(6,2,0);
530 localKeys_[13] = LocalKey(7,2,0);
531 localKeys_[14] = LocalKey(8,2,0);
532
533 // Quadrilateral sides shape functions
534 localKeys_[15] = LocalKey(0,1,0);
535 localKeys_[16] = LocalKey(1,1,0);
536 localKeys_[17] = LocalKey(2,1,0);
537
538 return;
539 }
540
541 // Now: the general case
542 DUNE_THROW(NotImplemented, "LagrangePrismLocalCoefficients not implemented for order " << k);
543
544 }
545
547 static constexpr std::size_t size ()
548 {
549 return binomial(k+2,2u) * (k+1);
550 }
551
553 const LocalKey& localKey (std::size_t i) const
554 {
555 return localKeys_[i];
556 }
557
558 private:
559 std::vector<LocalKey> localKeys_;
560 };
561
566 template<class LocalBasis>
567 class LagrangePrismLocalInterpolation
568 {
569 public:
570
578 template<typename F, typename C>
579 void interpolate (const F& ff, std::vector<C>& out) const
580 {
581 constexpr auto dim = LocalBasis::Traits::dimDomain;
582 constexpr auto k = LocalBasis::order();
583 using D = typename LocalBasis::Traits::DomainType;
584 using DF = typename LocalBasis::Traits::DomainFieldType;
585
586 auto&& f = Impl::makeFunctionWithCallOperator<D>(ff);
587
588 out.resize(LocalBasis::size());
589
590 // Specialization for zero-order case
591 if (k==0)
592 {
593 auto center = ReferenceElements<DF,dim>::general(GeometryTypes::prism).position(0,0);
594 out[0] = f(center);
595 return;
596 }
597
598 // Specialization for first-order case
599 if (k==1)
600 {
601 for (unsigned int i=0; i<LocalBasis::size(); i++)
602 {
603 auto vertex = ReferenceElements<DF,3>::general(GeometryTypes::prism).position(i,3);
604 out[i] = f(vertex);
605 }
606 return;
607 }
608
609 if (k==2)
610 {
611 out[0] = f( D( {0.0, 0.0, 0.0} ) );
612 out[1] = f( D( {1.0, 0.0, 0.0} ) );
613 out[2] = f( D( {0.0, 1.0, 0.0} ) );
614 out[3] = f( D( {0.0, 0.0, 1.0} ) );
615 out[4] = f( D( {1.0, 0.0, 1.0} ) );
616 out[5] = f( D( {0.0, 1.0, 1.0} ) );
617 out[6] = f( D( {0.0, 0.0, 0.5} ) );
618 out[7] = f( D( {1.0, 0.0, 0.5} ) );
619 out[8] = f( D( {0.0, 1.0, 0.5} ) );
620 out[9] = f( D( {0.5, 0.0, 0.0} ) );
621 out[10] = f( D( {0.0, 0.5, 0.0} ) );
622 out[11] = f( D( {0.5, 0.5, 0.0} ) );
623 out[12] = f( D( {0.5, 0.0, 1.0} ) );
624 out[13] = f( D( {0.0, 0.5, 1.0} ) );
625 out[14] = f( D( {0.5, 0.5, 1.0} ) );
626 out[15] = f( D( {0.5, 0.0, 0.5} ) );
627 out[16] = f( D( {0.0, 0.5, 0.5} ) );
628 out[17] = f( D( {0.5, 0.5, 0.5} ) );
629
630 return;
631 }
632
633 DUNE_THROW(NotImplemented, "LagrangePrismLocalInterpolation not implemented for order " << k);
634 }
635
636 };
637
638} } // namespace Dune::Impl
639
640namespace Dune
641{
648 template<class D, class R, int k>
649 class LagrangePrismLocalFiniteElement
650 {
651 public:
654 using Traits = LocalFiniteElementTraits<Impl::LagrangePrismLocalBasis<D,R,k>,
655 Impl::LagrangePrismLocalCoefficients<k>,
656 Impl::LagrangePrismLocalInterpolation<Impl::LagrangePrismLocalBasis<D,R,k> > >;
657
663 LagrangePrismLocalFiniteElement() {}
664
667 const typename Traits::LocalBasisType& localBasis () const
668 {
669 return basis_;
670 }
671
674 const typename Traits::LocalCoefficientsType& localCoefficients () const
675 {
676 return coefficients_;
677 }
678
681 const typename Traits::LocalInterpolationType& localInterpolation () const
682 {
683 return interpolation_;
684 }
685
687 static constexpr std::size_t size ()
688 {
689 return binomial(k+2,2) * (k+1);
690 }
691
694 static constexpr GeometryType type ()
695 {
696 return GeometryTypes::prism;
697 }
698
699 private:
700 Impl::LagrangePrismLocalBasis<D,R,k> basis_;
701 Impl::LagrangePrismLocalCoefficients<k> coefficients_;
702 Impl::LagrangePrismLocalInterpolation<Impl::LagrangePrismLocalBasis<D,R,k> > interpolation_;
703 };
704
705} // namespace Dune
706
707#endif // DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGEPRISM_HH
Dune::GeometryType
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:123
Dune::LagrangePrismLocalFiniteElement
Lagrange finite element for 3d prisms with arbitrary compile-time polynomial order.
Definition: lagrangeprism.hh:650
Dune::LagrangePrismLocalFiniteElement::size
static constexpr std::size_t size()
The number of shape functions.
Definition: lagrangeprism.hh:687
Dune::LagrangePrismLocalFiniteElement::localInterpolation
const Traits::LocalInterpolationType & localInterpolation() const
Returns object that evaluates degrees of freedom.
Definition: lagrangeprism.hh:681
Dune::LagrangePrismLocalFiniteElement::type
static constexpr GeometryType type()
The reference element that the local finite element is defined on.
Definition: lagrangeprism.hh:694
Dune::LagrangePrismLocalFiniteElement::localCoefficients
const Traits::LocalCoefficientsType & localCoefficients() const
Returns the assignment of the degrees of freedom to the element subentities.
Definition: lagrangeprism.hh:674
Dune::LagrangePrismLocalFiniteElement::LagrangePrismLocalFiniteElement
LagrangePrismLocalFiniteElement()
Default constructor.
Definition: lagrangeprism.hh:663
Dune::LagrangePrismLocalFiniteElement::localBasis
const Traits::LocalBasisType & localBasis() const
Returns the local basis, i.e., the set of shape functions.
Definition: lagrangeprism.hh:667
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::prism
constexpr GeometryType prism
GeometryType representing a 3D prism.
Definition: type.hh:540
Dune::GeometryTypes::vertex
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:504
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:289
math.hh
Some useful basic math stuff.
Dune
Dune namespace.
Definition: alignedallocator.hh:11
Dune::binomial
static constexpr T binomial(const T &n, const T &k) noexcept
calculate the binomial coefficient n over k as a constexpr
Definition: math.hh:128
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 (Dec 21, 23:30, 2024)