Dune Core Modules (2.9.0)

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