Dune Core Modules (2.7.1)

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 
20 namespace 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 
640 namespace 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:279
Dune::LagrangePrismLocalFiniteElement
Lagrange finite element for 3d prisms with arbitrary compile-time polynomial order.
Definition: lagrangeprism.hh:650
Dune::LagrangePrismLocalFiniteElement::localInterpolation
const Traits::LocalInterpolationType & localInterpolation() const
Returns object that evaluates degrees of freedom.
Definition: lagrangeprism.hh:681
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::size
static constexpr std::size_t size()
The number of shape functions.
Definition: lagrangeprism.hh:687
Dune::LagrangePrismLocalFiniteElement::localBasis
const Traits::LocalBasisType & localBasis() const
Returns the local basis, i.e., the set of shape functions.
Definition: lagrangeprism.hh:667
Dune::LagrangePrismLocalFiniteElement::type
static constexpr GeometryType type()
The reference element that the local finite element is defined on.
Definition: lagrangeprism.hh:694
Dune::LagrangePrismLocalFiniteElement::LagrangePrismLocalFiniteElement
LagrangePrismLocalFiniteElement()
Default constructor.
Definition: lagrangeprism.hh:663
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:833
Dune::GeometryTypes::vertex
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:797
Dune::Hybrid::accumulate
T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:290
math.hh
Some useful basic math stuff.
Dune
Dune namespace.
Definition: alignedallocator.hh:14
Dune::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: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.80.0 (May 16, 22:29, 2024)