Dune Core Modules (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 
21 namespace 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 
639 namespace 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::localInterpolation
const Traits::LocalInterpolationType & localInterpolation() const
Returns object that evaluates degrees of freedom.
Definition: lagrangeprism.hh:673
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::size
static constexpr std::size_t size()
The number of shape functions.
Definition: lagrangeprism.hh:679
Dune::LagrangePrismLocalFiniteElement::localBasis
const Traits::LocalBasisType & localBasis() const
Returns the local basis, i.e., the set of shape functions.
Definition: lagrangeprism.hh:659
Dune::LagrangePrismLocalFiniteElement::type
static constexpr GeometryType type()
The reference element that the local finite element is defined on.
Definition: lagrangeprism.hh:686
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: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
constexpr static 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.80.0 (May 2, 22:35, 2024)