DUNE-FEM (unstable)

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