Dune Core Modules (2.7.0)

lagrangepyramid.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_LAGRANGEPYRAMID_HH
4 #define DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGEPYRAMID_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 LagrangePyramidLocalBasis
33  {
34  public:
35  using Traits = LocalBasisTraits<D,3,FieldVector<D,3>,R,1,FieldVector<R,1>,FieldMatrix<R,1,3> >;
36 
39  static constexpr std::size_t size ()
40  {
41  std::size_t result = 0;
42  for (unsigned int i=0; i<=k; i++)
43  result += power(i+1,2);
44  return result;
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  if(in[0] > in[1])
63  {
64  out[0] = (1-in[0])*(1-in[1])-in[2]*(1-in[1]);
65  out[1] = in[0]*(1-in[1])-in[2]*in[1];
66  out[2] = (1-in[0])*in[1]-in[2]*in[1];
67  out[3] = in[0]*in[1]+in[2]*in[1];
68  }
69  else
70  {
71  out[0] = (1-in[0])*(1-in[1])-in[2]*(1-in[0]);
72  out[1] = in[0]*(1-in[1])-in[2]*in[0];
73  out[2] = (1-in[0])*in[1]-in[2]*in[0];
74  out[3] = in[0]*in[1]+in[2]*in[0];
75  }
76 
77  out[4] = in[2];
78 
79  return;
80  }
81 
82  if (k==2)
83  {
84  // transform to reference element with base [-1,1]^2
85  const R x = 2.0*in[0] + in[2] - 1.0;
86  const R y = 2.0*in[1] + in[2] - 1.0;
87  const R z = in[2];
88 
89  if (x > y)
90  {
91  // vertices
92  out[0] = 0.25*(x + z)*(x + z - 1)*(y - z - 1)*(y - z);
93  out[1] = -0.25*(x + z)*(y - z)*((x + z + 1)*(-y + z + 1) - 4*z) - z*(x - y);
94  out[2] = 0.25*(x + z)*(y - z)*(y - z + 1)*(x + z - 1);
95  out[3] = 0.25*(y - z)*(x + z)*(y - z + 1)*(x + z + 1);
96  out[4] = z*(2*z - 1);
97 
98  // lower edges
99  out[5] = -0.5*(y - z + 1)*(x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1));
100  out[6] = -0.5*(y - z + 1)*(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1));
101  out[7] = -0.5*(x + z - 1)*(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1));
102  out[8] = -0.5*(y - z + 1)*(x + z - 1)*(x + 1)*y;
103 
104  // upper edges
105  out[9] = z*(x + z - 1)*(y - z - 1);
106  out[10] = -z*((x + z + 1)*(y - z - 1) + 4*z);
107  out[11] = -z*(y - z + 1)*(x + z - 1);
108  out[12] = z*(y - z + 1)*(x + z + 1);
109 
110  // base face
111  out[13] = (y - z + 1)*(x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1));
112  }
113  else
114  {
115  // vertices
116  out[0] = 0.25*(y + z)*(y + z - 1)*(x - z - 1)*(x - z);
117  out[1] = -0.25*(x - z)*(y + z)*(x - z + 1)*(-y - z + 1);
118  out[2] = 0.25*(x - z)*(y + z)*((x - z - 1)*(y + z + 1) + 4*z) + z*(x - y);
119  out[3] = 0.25*(y + z)*(x - z)*(x - z + 1)*(y + z + 1);
120  out[4] = z*(2*z - 1);
121 
122  // lower edges
123  out[5] = -0.5*(y + z - 1)*(((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1));
124  out[6] = -0.5*(x - z + 1)*(y + z - 1)*(y + 1)*x;
125  out[7] = -0.5*(x - z + 1)*(y + z - 1)*(x - 1)*y;
126  out[8] = -0.5*(x - z + 1)*(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1));
127 
128  // upper edges
129  out[9] = z*(y + z - 1)*(x - z - 1);
130  out[10] = -z*(x - z + 1)*(y + z - 1);
131  out[11] = -z*((y + z + 1)*(x - z - 1) + 4*z);
132  out[12] = z*(x - z + 1)*(y + z + 1);
133 
134  // base face
135  out[13] = (x - z + 1)*(y + z - 1)*((y + 1)*(x - 1) - z*(x - y - z - 1));
136  }
137 
138  return;
139  }
140 
141  DUNE_THROW(NotImplemented, "LagrangePyramidLocalBasis::evaluateFunction for order " << k);
142  }
143 
149  void evaluateJacobian(const typename Traits::DomainType& in,
150  std::vector<typename Traits::JacobianType>& out) const
151  {
152  out.resize(size());
153 
154  // Specialization for k==0
155  if (k==0)
156  {
157  std::fill(out[0][0].begin(), out[0][0].end(), 0);
158  return;
159  }
160 
161  if (k==1)
162  {
163  if(in[0] > in[1])
164  {
165  out[0][0] = {-1 + in[1], -1 + in[0] + in[2], -1 + in[1]};
166  out[1][0] = { 1 - in[1], -in[0] - in[2], -in[1]};
167  out[2][0] = { -in[1], 1 - in[0] - in[2], -in[1]};
168  out[3][0] = { in[1], in[0] + in[2], in[1]};
169  }
170  else
171  {
172  out[0][0] = {-1 + in[1] + in[2], -1 + in[0], -1 + in[0]};
173  out[1][0] = { 1 - in[1] - in[2], -in[0], -in[0]};
174  out[2][0] = { -in[1] - in[2], 1 - in[0], -in[0]};
175  out[3][0] = { in[1] + in[2], in[0], in[0]};
176  }
177 
178  out[4][0] = {0, 0, 1};
179  return;
180  }
181 
182  if (k==2)
183  {
184  // transform to reference element with base [-1,1]^2
185  const R x = 2.0*in[0] + in[2] - 1.0;
186  const R y = 2.0*in[1] + in[2] - 1.0;
187  const R z = in[2];
188 
189  // transformation of the gradient leads to a multiplication
190  // with the Jacobian [2 0 0; 0 2 0; 1 1 1]
191  if (x > y)
192  {
193  // vertices
194  out[0][0][0] = 0.5*(y - z - 1)*(y - z)*(2*x + 2*z - 1);
195  out[0][0][1] = 0.5*(x + z)*(x + z - 1)*(2*y - 2*z - 1);
196  out[0][0][2] = 0.5*(out[0][0][0] + out[0][0][1])
197  + 0.25*((2*x + 2*z - 1)*(y - z - 1)*(y - z)
198  + (x + z)*(x + z - 1)*(-2*y + 2*z + 1));
199 
200  out[1][0][0] = 2*(-0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z)
201  + (x + z)*(y - z)*(-y + z + 1)) - z);
202  out[1][0][1] = 2*(-0.25*((x + z)*((x + z + 1)*(-y + z + 1) - 4*z)
203  + (x + z)*(y - z)*(-(x + z + 1))) + z);
204  out[1][0][2] = 0.5*(out[1][0][0] + out[1][0][1])
205  - 0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z)
206  - (x + z)*((x + z + 1)*(-y + z + 1) - 4*z)
207  + (x + z)*(y - z)*(x - y + 2*z - 2))
208  - (x - y);
209 
210  out[2][0][0] = 0.5*(y - z)*(y - z + 1)*(2*x + 2*z - 1);
211  out[2][0][1] = 0.5*(x + z)*(2*y - 2*z + 1)*(x + z - 1);
212  out[2][0][2] = 0.5*(out[2][0][0] + out[2][0][1])
213  + 0.25*((y - x - 2*z)*(y - z + 1)*(x + z - 1)
214  + (x + z)*(y - z)*(y - x - 2*z + 2));
215 
216  out[3][0][0] = 0.5*(y - z)*(2*x + 2*z + 1)*(y - z + 1);
217  out[3][0][1] = 0.5*(2*y - 2*z + 1)*(x + z)*(x + z + 1);
218  out[3][0][2] = 0.5*(out[3][0][0] + out[3][0][1])
219  + 0.25*((y - x - 2*z)*(y - z + 1)*(x + z + 1)
220  + (y - z)*(x + z)*(y - x - 2*z));
221 
222  out[4][0][0] = 0;
223  out[4][0][1] = 0;
224  out[4][0][2] = 4*z - 1;
225 
226  // lower edges
227  out[5][0][0] = -((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
228  + (y - z + 1)*(x + z - 1)*((y - 1) + z));
229  out[5][0][1] = -((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
230  + (y - z + 1)*(x + z - 1)*((x + 1) - z));
231  out[5][0][2] = 0.5*(out[5][0][0] + out[5][0][1])
232  - 0.5*((-x + y - 2*z + 2)*((y - 1)*(x + 1) + z*(x - y + z + 1))
233  + (y - z + 1)*(x + z - 1)*(x - y + 2*z + 1));
234 
235  out[6][0][0] = -(y - z + 1)*(2*x + z + 1)*(y - 1);
236  out[6][0][1] = -(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1)
237  + (y - z + 1)*((x + z + 1)*x + 2*z));
238  out[6][0][2] = 0.5*(out[6][0][0] + out[6][0][1])
239  - 0.5*(-(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1))
240  + (y - z + 1)*(((y - 1)*x - 1) + 2*y + 1));
241 
242  out[7][0][0] = -(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1)
243  + (x + z - 1)*((y - z - 1)*y + 2*z));
244  out[7][0][1] = -(x + z - 1)*(2*y - z - 1)*(x + 1);
245  out[7][0][2] = 0.5*(out[7][0][0] + out[7][0][1])
246  - 0.5*(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1)
247  + (x + z - 1)*((-(x + 1)*y - 1) + 2*x + 1));
248 
249  out[8][0][0] = -(y - z + 1)*(2*x + z)*y;
250  out[8][0][1] = -(2*y - z + 1)*(x + z - 1)*(x + 1);
251  out[8][0][2] = 0.5*(out[8][0][0] + out[8][0][1])
252  - 0.5*(-x + y - 2*z + 2)*(x + 1)*y;
253 
254  // upper edges
255  out[9][0][0] = 2*z*(y - z - 1);
256  out[9][0][1] = 2*z*(x + z - 1);
257  out[9][0][2] = 0.5*(out[9][0][0] + out[9][0][1])
258  + (x + z - 1)*(y - z - 1) + z*(-x + y - 2*z);
259 
260  out[10][0][0] = -2*z*(y - z - 1);
261  out[10][0][1] = -2*z*(x + z + 1);
262  out[10][0][2] = 0.5*(out[10][0][0] + out[10][0][1])
263  - ((x + z + 1)*(y - z - 1) + 4*z)
264  - z*(-x + y - 2*z + 2);
265 
266  out[11][0][0] = -2*z*(y - z + 1);
267  out[11][0][1] = -2*z*(x + z - 1);
268  out[11][0][2] = 0.5*(out[11][0][0] + out[11][0][1])
269  - (y - z + 1)*(x + z - 1) - z*(-x + y - 2*z + 2);
270 
271  out[12][0][0] = 2*z*(y - z + 1);
272  out[12][0][1] = 2*z*(x + z + 1);
273  out[12][0][2] = 0.5*(out[12][0][0] + out[12][0][1])
274  + (y - z + 1)*(x + z + 1) + z*(-x + y - 2*z);
275 
276  // base face
277  out[13][0][0] = 2*((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
278  + (y - z + 1)*(x + z - 1)*(y - 1 + z));
279  out[13][0][1] = 2*((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
280  + (y - z + 1)*(x + z - 1)*(x + 1 - z));
281  out[13][0][2] = 0.5*(out[13][0][0] + out[13][0][1])
282  + ((-x + y - 2*z + 2)*((y - 1)*(x + 1) + z*(x - y + z + 1))
283  + (y - z + 1)*(x + z - 1)*(x - y + 2*z + 1));
284  }
285  else
286  {
287  // vertices
288  out[0][0][0] = 0.5*(y + z)*(y + z - 1)*(2*x - 2*z - 1);
289  out[0][0][1] = 0.5*(2*y + 2*z - 1)*(x - z - 1)*(x - z);
290  out[0][0][2] = 0.5*(out[0][0][0] + out[0][0][1])
291  + 0.25*((2*y + 2*z - 1)*(x - z - 1)*(x - z)
292  + (y + z)*(y + z - 1)*(-2*x + 2*z + 1));
293 
294  out[1][0][0] = -0.5*(y + z)*(2*x - 2*z + 1)*(-y - z + 1);
295  out[1][0][1] = -0.5*(x - z)*(x - z + 1)*(-2*y - 2*z + 1);
296  out[1][0][2] = 0.5*(out[1][0][0] + out[1][0][1])
297  - 0.25*((x - y - 2*z)*(x - z + 1)*(-y - z + 1)
298  + (x - z)*(y + z)*(-x + y + 2*z - 2));
299 
300  out[2][0][0] = 0.5*((y + z)*((x - z - 1)*(y + z + 1) + 4*z)
301  + (x - z)*(y + z)*(y + z + 1) + 4*z);
302  out[2][0][1] = 0.5*((x - z)*((x - z - 1)*(y + z + 1) + 4*z)
303  + (x - z)*(y + z)*(x - z - 1) - 4*z);
304  out[2][0][2] = 0.5*(out[2][0][0] + out[2][0][1])
305  + 0.25*((x - y - 2*z)*((x - z - 1)*(y + z + 1) + 4*z)
306  + (x - z)*(y + z)*(x - y - 2*z + 2) + 4*(x - y));
307 
308  out[3][0][0] = 0.5*(y + z)*(2*x - 2*z + 1)*(y + z + 1);
309  out[3][0][1] = 0.5*(x - z)*(x - z + 1)*(2*y + 2*z + 1);
310  out[3][0][2] = 0.5*(out[3][0][0] + out[3][0][1])
311  + 0.25*((x - y - 2*z)*(x - z + 1)*(y + z + 1)
312  + (y + z)*(x - z)*(x - y - 2*z));
313 
314  out[4][0][0] = 0;
315  out[4][0][1] = 0;
316  out[4][0][2] = 4*z - 1;
317 
318  // lower edges
319  out[5][0][0] = -(y + z - 1)*(2*x - z - 1)*(y + 1);
320  out[5][0][1] = -(((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1)
321  + (y + z - 1)*((x - z - 1)*x + 2*z));
322  out[5][0][2] = 0.5*(out[5][0][0] + out[5][0][1])
323  - 0.5*((((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1))
324  + (y + z - 1)*((-(y + 1)*x - 1) + 2*y + 1));
325 
326  out[6][0][0] = -(2*x - z + 1)*(y + z - 1)*(y + 1);
327  out[6][0][1] = -(x - z + 1)*(2*y + z)*x;
328  out[6][0][2] = 0.5*(out[6][0][0] + out[6][0][1])
329  - 0.5*(x - y - 2*z + 2)*(y + 1)*x;
330 
331  out[7][0][0] = -(2*x - z)*(y + z - 1)*y;
332  out[7][0][1] = -(x - z + 1)*(2*y + z - 1)*(x - 1);
333  out[7][0][2] = 0.5*(out[7][0][0] + out[7][0][1])
334  - 0.5*(x - y - 2*z + 2)*(x - 1)*y;
335 
336  out[8][0][0] = -(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1)
337  + (x - z + 1)*((y + z + 1)*y + 2*z));
338  out[8][0][1] = -(x - z + 1)*(2*y + z + 1)*(x - 1);
339  out[8][0][2] = 0.5*(out[8][0][0] + out[8][0][1])
340  - 0.5*(-(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1))
341  + (x - z + 1)*(((x - 1)*y - 1) + 2*x + 1));
342 
343  // upper edges
344  out[9][0][0] = 2*z*(y + z - 1);
345  out[9][0][1] = 2*z*(x - z - 1);
346  out[9][0][2] = 0.5*(out[9][0][0] + out[9][0][1])
347  + (y + z - 1)*(x - z - 1) + z*(x - y - 2*z);
348 
349  out[10][0][0] = -2*z*(y + z - 1);
350  out[10][0][1] = -2*z*(x - z + 1);
351  out[10][0][2] = 0.5*(out[10][0][0] + out[10][0][1])
352  - (x - z + 1)*(y + z - 1) - z*(x - y - 2*z + 2);
353 
354  out[11][0][0] = -2*z*(y + z + 1);
355  out[11][0][1] = -2*z*(x - z - 1);
356  out[11][0][2] = 0.5*(out[11][0][0] + out[11][0][1])
357  - ((y + z + 1)*(x - z - 1) + 4*z) - z*(x - y - 2*z + 2);
358 
359  out[12][0][0] = 2*z*(y + z + 1);
360  out[12][0][1] = 2*z*(x - z + 1);
361  out[12][0][2] = 0.5*(out[12][0][0] + out[12][0][1])
362  + (x - z + 1)*(y + z + 1) + z*(x - y - 2*z);
363 
364  // base face
365  out[13][0][0] = 2*((y + z - 1)*((y + 1)*(x - 1) - z*(x - y - z - 1))
366  + (x - z + 1)*(y + z - 1)*(y + 1 - z));
367  out[13][0][1] = 2*((x - z + 1)*((y + 1)*(x - 1) - z*(x - y - z - 1))
368  + (x - z + 1)*(y + z - 1)*(x - 1 + z));
369  out[13][0][2] = 0.5*(out[13][0][0] + out[13][0][1])
370  + (x - y - 2*z + 2)*((y + 1)*(x - 1) - z*(x - y - z - 1))
371  + (x - z + 1)*(y + z - 1)*(-(x - y - 2*z - 1));
372  }
373 
374  return;
375  }
376 
377  DUNE_THROW(NotImplemented, "LagrangePyramidLocalBasis::evaluateJacobian for order " << k);
378  }
379 
386  void partial(const std::array<unsigned int,3>& order,
387  const typename Traits::DomainType& in,
388  std::vector<typename Traits::RangeType>& out) const
389  {
390  auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
391 
392  out.resize(size());
393 
394  if (totalOrder == 0)
395  {
396  evaluateFunction(in, out);
397  return;
398  }
399 
400  if (k==0)
401  {
402  out[0] = 0;
403  return;
404  }
405 
406  if (k==1)
407  {
408  if (totalOrder == 1)
409  {
410  auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
411  if (in[0] > in[1])
412  {
413  switch (direction)
414  {
415  case 0:
416  out[0] = -1 + in[1];
417  out[1] = 1 - in[1];
418  out[2] = -in[1];
419  out[3] = in[1];
420  out[4] = 0;
421  break;
422  case 1:
423  out[0] = -1 + in[0] + in[2];
424  out[1] = -in[0] - in[2];
425  out[2] = 1 - in[0] - in[2];
426  out[3] = in[0]+in[2];
427  out[4] = 0;
428  break;
429  case 2:
430  out[0] = -1 + in[1];
431  out[1] = -in[1];
432  out[2] = -in[1];
433  out[3] = in[1];
434  out[4] = 1;
435  break;
436  default:
437  DUNE_THROW(RangeError, "Component out of range.");
438  }
439  }
440  else /* (in[0] <= in[1]) */
441  {
442  switch (direction)
443  {
444  case 0:
445  out[0] = -1 + in[1] + in[2];
446  out[1] = 1 - in[1] - in[2];
447  out[2] = -in[1] - in[2];
448  out[3] = in[1] + in[2];
449  out[4] = 0;
450  break;
451  case 1:
452  out[0] = -1 + in[0];
453  out[1] = -in[0];
454  out[2] = 1 - in[0];
455  out[3] = in[0];
456  out[4] = 0;
457  break;
458  case 2:
459  out[0] = -1 + in[0];
460  out[1] = -in[0];
461  out[2] = -in[0];
462  out[3] = in[0];
463  out[4] = 1;
464  break;
465  default:
466  DUNE_THROW(RangeError, "Component out of range.");
467  }
468  }
469  } else if (totalOrder == 2)
470  {
471  if ((order[0] == 1 && order[1] == 1) ||
472  (order[1] == 1 && order[2] == 1 && in[0] > in[1]) ||
473  (order[0] == 1 && order[2] == 1 && in[0] <=in[1]))
474  {
475  out = {1, -1, -1, 1, 0};
476  } else
477  {
478  out = {0, 0, 0, 0, 0};
479  }
480 
481  } else
482  {
483  out = {0, 0, 0, 0, 0};
484  }
485 
486  return;
487  }
488 
489  if (k==2)
490  {
491  if (totalOrder == 1)
492  {
493  // transform to reference element with base [-1,1]^2
494  const R x = 2.0*in[0] + in[2] - 1.0;
495  const R y = 2.0*in[1] + in[2] - 1.0;
496  const R z = in[2];
497 
498  auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
499 
500  // transformation of the gradient leads to a multiplication
501  // with the Jacobian [2 0 0; 0 2 0; 1 1 1]
502  if (x > y)
503  {
504  switch (direction)
505  {
506  case 0:
507  out[0] = 0.5*(y - z - 1)*(y - z)*(2*x + 2*z - 1);
508  out[1] = 2*(-0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z) + (x + z)*(y - z)*(-y + z + 1)) - z);
509  out[2] = 0.5*(y - z)*(y - z + 1)*(2*x + 2*z - 1);
510  out[3] = 0.5*(y - z)*(2*x + 2*z + 1)*(y - z + 1);
511  out[4] = 0;
512  out[5] = -((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*((y - 1) + z));
513  out[6] = -(y - z + 1)*(2*x + z + 1)*(y - 1);
514  out[7] = -(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1) + (x + z - 1)*((y - z - 1)*y + 2*z));
515  out[8] = -(y - z + 1)*(2*x + z)*y;
516  out[9] = 2*z*(y - z - 1);
517  out[10] = -2*z*(y - z - 1);
518  out[11] = -2*z*(y - z + 1);
519  out[12] = 2*z*(y - z + 1);
520  out[13] = 2*((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*(y - 1 + z));
521  break;
522  case 1:
523  out[0] = 0.5*(x + z)*(x + z - 1)*(2*y - 2*z - 1);
524  out[1] = 2*(-0.25*((x + z)*((x + z + 1)*(-y + z + 1) - 4*z) + (x + z)*(y - z)*(-(x + z + 1))) + z);
525  out[2] = 0.5*(x + z)*(2*y - 2*z + 1)*(x + z - 1);
526  out[3] = 0.5*(2*y - 2*z + 1)*(x + z)*(x + z + 1);
527  out[4] = 0;
528  out[5] = -((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*((x + 1) - z));
529  out[6] = -(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1) + (y - z + 1)*((x + z + 1)*x + 2*z));
530  out[7] = -(x + z - 1)*(2*y - z - 1)*(x + 1);
531  out[8] = -(2*y - z + 1)*(x + z - 1)*(x + 1);
532  out[9] = 2*z*(x + z - 1);
533  out[10] = -2*z*(x + z + 1);
534  out[11] = -2*z*(x + z - 1);
535  out[12] = 2*z*(x + z + 1);
536  out[13] = 2*((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*(x + 1 - z));
537  break;
538  case 2:
539  out[0] = -((y - z)*(2*x + 2*z - 1)*(z - y + 1))/2;
540  out[1] = ((y - z + 1)*(y - 2*x + z + 2*x*y - 2*x*z + 2*y*z - 2*z*z))/2;
541  out[2] = ((y - z)*(2*x + 2*z - 1)*(y - z + 1))/2;
542  out[3] = ((y - z)*(2*x + 2*z + 1)*(y - z + 1))/2;
543  out[4] = 4*z - 1;
544  out[5] = -((y - z + 1)*(2*y - 3*x + z + 2*x*y + 6*x*z - 2*y*z + 2*x*x + 4*z*z - 3))/2;
545  out[6] = -((y - z + 1)*(3*y - 2*x + z + 3*x*y + x*z + y*z + x*x - 1))/2;
546  out[7] = z - z*(2*x + 1) - ((2*z - y*(z - y + 1))*(x + z - 1))/2 - ((2*x - y*(x + 1))*(x + z - 1))/2 + ((x + 1)*(x + z - 1)*(z - 2*y + 1))/2 + y*(x + 1)*(z - y + 1);
547  out[8] = -((y - z + 1)*(y + z + 3*x*y + x*z + y*z + x*x - 1))/2;
548  out[9] = -(x + 3*z - 1)*(z - y + 1);
549  out[10] = (x + z + 1)*(z - y + 1) - 2*y*z - 6*z + 2*z*z;
550  out[11] = -(x + 3*z - 1)*(y - z + 1);
551  out[12] = (x + 3*z + 1)*(y - z + 1);
552  out[13] = (y - z + 1)*(2*y - 3*x + z + 2*x*y + 6*x*z - 2*y*z + 2*x*x + 4*z*z - 3);
553  break;
554  default:
555  DUNE_THROW(RangeError, "Component out of range.");
556  }
557  }
558  else // x <= y
559  {
560  switch (direction)
561  {
562  case 0:
563  out[0] = -((y + z)*(2*z - 2*x + 1)*(y + z - 1))/2;
564  out[1] = ((y + z)*(2*x - 2*z + 1)*(y + z - 1))/2;
565  out[2] = -((y + z + 1)*(y - 3*z - 2*x*y - 2*x*z + 2*y*z + 2*z*z))/2;
566  out[3] = ((y + z)*(2*x - 2*z + 1)*(y + z + 1))/2;
567  out[4] = 0;
568  out[5] = (y + 1)*(y + z - 1)*(z - 2*x + 1);
569  out[6] = -(y + 1)*(2*x - z + 1)*(y + z - 1);
570  out[7] = -y*(2*x - z)*(y + z - 1);
571  out[8] = z - z*(2*x + 1) - (2*z + y*(y + z + 1))*(x - z + 1) - y*(x - 1)*(y + z + 1);
572  out[9] = 2*z*(y + z - 1);
573  out[10] = -2*z*(y + z - 1);
574  out[11] = -2*z*(y + z + 1);
575  out[12] = 2*z*(y + z + 1);
576  out[13] = 2*(y + z - 1)*(2*x - z + 2*x*y - 2*x*z + 2*z*z);
577  break;
578  case 1:
579  out[0] = -(x - z)*(y + z - 0.5)*(z - x + 1);
580  out[1] = ((x - z)*(2*y + 2*z - 1)*(x - z + 1))/2;
581  out[2] = -((z - x + 1)*(x + 3*z + 2*x*y + 2*x*z - 2*y*z - 2*z*z))/2;
582  out[3] = ((x - z)*(2*y + 2*z + 1)*(x - z + 1))/2;
583  out[4] = 0;
584  out[5] = z - z*(2*y + 1) - (2*z - x*(z - x + 1))*(y + z - 1) + x*(y + 1)*(z - x + 1);
585  out[6] = -x*(2*y + z)*(x - z + 1);
586  out[7] = -(x - 1)*(x - z + 1)*(2*y + z - 1);
587  out[8] = -(x - 1)*(x - z + 1)*(2*y + z + 1);
588  out[9] = -2*z*(z - x + 1);
589  out[10] = -2*z*(x - z + 1);
590  out[11] = 2*z*(z - x + 1);
591  out[12] = 2*z*(x - z + 1);
592  out[13] = 2*(x - z + 1)*(2*x*y - z - 2*y + 2*y*z + 2*z*z);
593  break;
594  case 2:
595  out[0] = -((x - z)*(2*y + 2*z - 1)*(z - x + 1))/2;
596  out[1] = ((x - z)*(2*y + 2*z - 1)*(x - z + 1))/2;
597  out[2] = ((x - z + 1)*(x - 2*y + z + 2*x*y + 2*x*z - 2*y*z - 2*z*z))/2;
598  out[3] = ((x - z)*(2*y + 2*z + 1)*(x - z + 1))/2;
599  out[4] = 4*z - 1;
600  out[5] = z - z*(2*y + 1) - ((2*z - x*(z - x + 1))*(y + z - 1))/2 - ((2*y - x*(y + 1))*(y + z - 1))/2 + ((y + 1)*(y + z - 1)*(z - 2*x + 1))/2 + x*(y + 1)*(z - x + 1);
601  out[6] = -((x - z + 1)*(x + z + 3*x*y + x*z + y*z + y*y - 1))/2;
602  out[7] = -((x - z + 1)*(3*x*y - 4*y - z - x + x*z + y*z + y*y + 1))/2;
603  out[8] = -((x - z + 1)*(3*x - 2*y + z + 3*x*y + x*z + y*z + y*y - 1))/2;
604  out[9] = -(z - x + 1)*(y + 3*z - 1);
605  out[10] = -(x - z + 1)*(y + 3*z - 1);
606  out[11] = (y + z + 1)*(z - x + 1) - 2*x*z - 6*z + 2*z*z;
607  out[12] = (x - z + 1)*(y + 3*z + 1);
608  out[13] = (x - z + 1)*(2*x - 3*y + z + 2*x*y - 2*x*z + 6*y*z + 2*y*y + 4*z*z - 3);
609  break;
610  default:
611  DUNE_THROW(RangeError, "Component out of range.");
612  }
613  }
614  } else {
615  DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
616  }
617 
618  return;
619  }
620 
621  DUNE_THROW(NotImplemented, "LagrangePyramidLocalBasis::partial for order " << k);
622  }
623 
625  static constexpr unsigned int order ()
626  {
627  return k;
628  }
629  };
630 
635  template<unsigned int k>
636  class LagrangePyramidLocalCoefficients
637  {
638  public:
640  LagrangePyramidLocalCoefficients ()
641  : localKeys_(size())
642  {
643  if (k==0)
644  {
645  localKeys_[0] = LocalKey(0,0,0);
646  return;
647  }
648 
649  if (k==1)
650  {
651  for (std::size_t i=0; i<size(); i++)
652  localKeys_[i] = LocalKey(i,3,0);
653  return;
654  }
655 
656  if (k==2)
657  {
658  // Vertex shape functions
659  localKeys_[0] = LocalKey(0,3,0);
660  localKeys_[1] = LocalKey(1,3,0);
661  localKeys_[2] = LocalKey(2,3,0);
662  localKeys_[3] = LocalKey(3,3,0);
663  localKeys_[4] = LocalKey(4,3,0);
664 
665  // Edge shape functions
666  localKeys_[5] = LocalKey(0,2,0);
667  localKeys_[6] = LocalKey(1,2,0);
668  localKeys_[7] = LocalKey(2,2,0);
669  localKeys_[8] = LocalKey(3,2,0);
670  localKeys_[9] = LocalKey(4,2,0);
671  localKeys_[10] = LocalKey(5,2,0);
672  localKeys_[11] = LocalKey(6,2,0);
673  localKeys_[12] = LocalKey(7,2,0);
674 
675  // base face shape function
676  localKeys_[13] = LocalKey(0,1,0);
677 
678  return;
679  }
680 
681  // No general case
682  DUNE_THROW(NotImplemented, "LagrangePyramidLocalCoefficients for order " << k);
683 
684  }
685 
687  static constexpr std::size_t size ()
688  {
689  std::size_t result = 0;
690  for (unsigned int i=0; i<=k; i++)
691  result += power(i+1,2);
692  return result;
693  }
694 
696  const LocalKey& localKey (std::size_t i) const
697  {
698  return localKeys_[i];
699  }
700 
701  private:
702  std::vector<LocalKey> localKeys_;
703  };
704 
709  template<class LocalBasis>
710  class LagrangePyramidLocalInterpolation
711  {
712  public:
713 
721  template<typename F, typename C>
722  void interpolate (const F& ff, std::vector<C>& out) const
723  {
724  constexpr auto k = LocalBasis::order();
725  using D = typename LocalBasis::Traits::DomainType;
726  using DF = typename LocalBasis::Traits::DomainFieldType;
727 
728  auto&& f = Impl::makeFunctionWithCallOperator<D>(ff);
729 
730  out.resize(LocalBasis::size());
731 
732  // Specialization for zero-order case
733  if (k==0)
734  {
735  auto center = ReferenceElements<DF,3>::general(GeometryTypes::pyramid).position(0,0);
736  out[0] = f(center);
737  return;
738  }
739 
740  // Specialization for first-order case
741  if (k==1)
742  {
743  for (unsigned int i=0; i<LocalBasis::size(); i++)
744  {
745  auto vertex = ReferenceElements<DF,3>::general(GeometryTypes::pyramid).position(i,3);
746  out[i] = f(vertex);
747  }
748  return;
749  }
750 
751  // Specialization for second-order case
752  if (k==2)
753  {
754  out[0] = f( D( {0.0, 0.0, 0.0} ) );
755  out[1] = f( D( {1.0, 0.0, 0.0} ) );
756  out[2] = f( D( {0.0, 1.0, 0.0} ) );
757  out[3] = f( D( {1.0, 1.0, 0.0} ) );
758  out[4] = f( D( {0.0, 0.0, 1.0} ) );
759  out[5] = f( D( {0.0, 0.5, 0.0} ) );
760  out[6] = f( D( {1.0, 0.5, 0.0} ) );
761  out[7] = f( D( {0.5, 0.0, 0.0} ) );
762  out[8] = f( D( {0.5, 1.0, 0.0} ) );
763  out[9] = f( D( {0.0, 0.0, 0.5} ) );
764  out[10] = f( D( {0.5, 0.0, 0.5} ) );
765  out[11] = f( D( {0.0, 0.5, 0.5} ) );
766  out[12] = f( D( {0.5, 0.5, 0.5} ) );
767  out[13] = f( D( {0.5, 0.5, 0.0} ) );
768 
769  return;
770  }
771 
772  DUNE_THROW(NotImplemented, "LagrangePyramidLocalInterpolation not implemented for order " << k);
773  }
774 
775  };
776 
777 } } // namespace Dune::Impl
778 
779 namespace Dune
780 {
787  template<class D, class R, int k>
788  class LagrangePyramidLocalFiniteElement
789  {
790  public:
793  using Traits = LocalFiniteElementTraits<Impl::LagrangePyramidLocalBasis<D,R,k>,
794  Impl::LagrangePyramidLocalCoefficients<k>,
795  Impl::LagrangePyramidLocalInterpolation<Impl::LagrangePyramidLocalBasis<D,R,k> > >;
796 
802  LagrangePyramidLocalFiniteElement() {}
803 
806  const typename Traits::LocalBasisType& localBasis () const
807  {
808  return basis_;
809  }
810 
813  const typename Traits::LocalCoefficientsType& localCoefficients () const
814  {
815  return coefficients_;
816  }
817 
820  const typename Traits::LocalInterpolationType& localInterpolation () const
821  {
822  return interpolation_;
823  }
824 
826  static constexpr std::size_t size ()
827  {
828  return Impl::LagrangePyramidLocalBasis<D,R,k>::size();
829  }
830 
833  static constexpr GeometryType type ()
834  {
835  return GeometryTypes::pyramid;
836  }
837 
838  private:
839  Impl::LagrangePyramidLocalBasis<D,R,k> basis_;
840  Impl::LagrangePyramidLocalCoefficients<k> coefficients_;
841  Impl::LagrangePyramidLocalInterpolation<Impl::LagrangePyramidLocalBasis<D,R,k> > interpolation_;
842  };
843 
844 } // namespace Dune
845 
846 #endif // DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGEPYRAMID_HH
Dune::GeometryType
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:279
Dune::LagrangePyramidLocalFiniteElement
Lagrange finite element for 3d pyramids with arbitrary compile-time polynomial order.
Definition: lagrangepyramid.hh:789
Dune::LagrangePyramidLocalFiniteElement::localCoefficients
const Traits::LocalCoefficientsType & localCoefficients() const
Returns the assignment of the degrees of freedom to the element subentities.
Definition: lagrangepyramid.hh:813
Dune::LagrangePyramidLocalFiniteElement::size
static constexpr std::size_t size()
The number of shape functions.
Definition: lagrangepyramid.hh:826
Dune::LagrangePyramidLocalFiniteElement::LagrangePyramidLocalFiniteElement
LagrangePyramidLocalFiniteElement()
Default constructor.
Definition: lagrangepyramid.hh:802
Dune::LagrangePyramidLocalFiniteElement::type
static constexpr GeometryType type()
The reference element that the local finite element is defined on.
Definition: lagrangepyramid.hh:833
Dune::LagrangePyramidLocalFiniteElement::localBasis
const Traits::LocalBasisType & localBasis() const
Returns the local basis, i.e., the set of shape functions.
Definition: lagrangepyramid.hh:806
Dune::LagrangePyramidLocalFiniteElement::localInterpolation
const Traits::LocalInterpolationType & localInterpolation() const
Returns object that evaluates degrees of freedom.
Definition: lagrangepyramid.hh:820
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::pyramid
constexpr GeometryType pyramid
GeometryType representing a 3D pyramid.
Definition: type.hh:827
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::power
constexpr Mantissa power(Mantissa m, Exponent p)
Power method for integer exponents.
Definition: math.hh:73
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)