Dune Core Modules (2.7.1)

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;
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)*(2*x + z - 1);
228  out[5][0][1] = -(x + z - 1)*(y - 1)*x - (y - z + 1)*(x + z - 1)*x;
229  out[5][0][2] = 0.5*(out[5][0][0] + out[5][0][1])
230  + 0.5*(x + z - 1)*(y - 1)*x - 0.5*(y - z + 1)*(y - 1)*x;
231 
232  out[6][0][0] = -(y - z + 1)*(2*x + z + 1)*(y - 1);
233  out[6][0][1] = -(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1)
234  + (y - z + 1)*((x + z + 1)*x + 2*z));
235  out[6][0][2] = 0.5*(out[6][0][0] + out[6][0][1])
236  - 0.5*(-(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1))
237  + (y - z + 1)*(((y - 1)*x - 1) + 2*y + 1));
238 
239  out[7][0][0] = -(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1)
240  + (x + z - 1)*((y - z - 1)*y + 2*z));
241  out[7][0][1] = -(x + z - 1)*(2*y - z - 1)*(x + 1);
242  out[7][0][2] = 0.5*(out[7][0][0] + out[7][0][1])
243  - 0.5*(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1)
244  + (x + z - 1)*((-(x + 1)*y - 1) + 2*x + 1));
245 
246  out[8][0][0] = -(y - z + 1)*(2*x + z)*y;
247  out[8][0][1] = -(2*y - z + 1)*(x + z - 1)*(x + 1);
248  out[8][0][2] = 0.5*(out[8][0][0] + out[8][0][1])
249  - 0.5*(-x + y - 2*z + 2)*(x + 1)*y;
250 
251  // upper edges
252  out[9][0][0] = 2*z*(y - z - 1);
253  out[9][0][1] = 2*z*(x + z - 1);
254  out[9][0][2] = 0.5*(out[9][0][0] + out[9][0][1])
255  + (x + z - 1)*(y - z - 1) + z*(-x + y - 2*z);
256 
257  out[10][0][0] = -2*z*(y - z - 1);
258  out[10][0][1] = -2*z*(x + z + 1);
259  out[10][0][2] = 0.5*(out[10][0][0] + out[10][0][1])
260  - ((x + z + 1)*(y - z - 1) + 4*z)
261  - z*(-x + y - 2*z + 2);
262 
263  out[11][0][0] = -2*z*(y - z + 1);
264  out[11][0][1] = -2*z*(x + z - 1);
265  out[11][0][2] = 0.5*(out[11][0][0] + out[11][0][1])
266  - (y - z + 1)*(x + z - 1) - z*(-x + y - 2*z + 2);
267 
268  out[12][0][0] = 2*z*(y - z + 1);
269  out[12][0][1] = 2*z*(x + z + 1);
270  out[12][0][2] = 0.5*(out[12][0][0] + out[12][0][1])
271  + (y - z + 1)*(x + z + 1) + z*(-x + y - 2*z);
272 
273  // base face
274  out[13][0][0] = 2*((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
275  + (y - z + 1)*(x + z - 1)*(y - 1 + z));
276  out[13][0][1] = 2*((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
277  + (y - z + 1)*(x + z - 1)*(x + 1 - z));
278  out[13][0][2] = 0.5*(out[13][0][0] + out[13][0][1])
279  + ((-x + y - 2*z + 2)*((y - 1)*(x + 1) + z*(x - y + z + 1))
280  + (y - z + 1)*(x + z - 1)*(x - y + 2*z + 1));
281  }
282  else
283  {
284  // vertices
285  out[0][0][0] = 0.5*(y + z)*(y + z - 1)*(2*x - 2*z - 1);
286  out[0][0][1] = 0.5*(2*y + 2*z - 1)*(x - z - 1)*(x - z);
287  out[0][0][2] = 0.5*(out[0][0][0] + out[0][0][1])
288  + 0.25*((2*y + 2*z - 1)*(x - z - 1)*(x - z)
289  + (y + z)*(y + z - 1)*(-2*x + 2*z + 1));
290 
291  out[1][0][0] = -0.5*(y + z)*(2*x - 2*z + 1)*(-y - z + 1);
292  out[1][0][1] = -0.5*(x - z)*(x - z + 1)*(-2*y - 2*z + 1);
293  out[1][0][2] = 0.5*(out[1][0][0] + out[1][0][1])
294  - 0.25*((x - y - 2*z)*(x - z + 1)*(-y - z + 1)
295  + (x - z)*(y + z)*(-x + y + 2*z - 2));
296 
297  out[2][0][0] = 0.5*((y + z)*((x - z - 1)*(y + z + 1) + 4*z)
298  + (x - z)*(y + z)*(y + z + 1) + 4*z);
299  out[2][0][1] = 0.5*((x - z)*((x - z - 1)*(y + z + 1) + 4*z)
300  + (x - z)*(y + z)*(x - z - 1) - 4*z);
301  out[2][0][2] = 0.5*(out[2][0][0] + out[2][0][1])
302  + 0.25*((x - y - 2*z)*((x - z - 1)*(y + z + 1) + 4*z)
303  + (x - z)*(y + z)*(x - y - 2*z + 2) + 4*(x - y));
304 
305  out[3][0][0] = 0.5*(y + z)*(2*x - 2*z + 1)*(y + z + 1);
306  out[3][0][1] = 0.5*(x - z)*(x - z + 1)*(2*y + 2*z + 1);
307  out[3][0][2] = 0.5*(out[3][0][0] + out[3][0][1])
308  + 0.25*((x - y - 2*z)*(x - z + 1)*(y + z + 1)
309  + (y + z)*(x - z)*(x - y - 2*z));
310 
311  out[4][0][0] = 0;
312  out[4][0][1] = 0;
313  out[4][0][2] = 4*z - 1;
314 
315  // lower edges
316  out[5][0][0] = -(y + z - 1)*(2*x - z - 1)*(y + 1);
317  out[5][0][1] = -(((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1)
318  + (y + z - 1)*((x - z - 1)*x + 2*z));
319  out[5][0][2] = 0.5*(out[5][0][0] + out[5][0][1])
320  - 0.5*((((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1))
321  + (y + z - 1)*((-(y + 1)*x - 1) + 2*y + 1));
322 
323  out[6][0][0] = -(2*x - z + 1)*(y + z - 1)*(y + 1);
324  out[6][0][1] = -(x - z + 1)*(2*y + z)*x;
325  out[6][0][2] = 0.5*(out[6][0][0] + out[6][0][1])
326  - 0.5*(x - y - 2*z + 2)*(y + 1)*x;
327 
328  out[7][0][0] = -(2*x - z)*(y + z - 1)*y;
329  out[7][0][1] = -(x - z + 1)*(2*y + z - 1)*(x - 1);
330  out[7][0][2] = 0.5*(out[7][0][0] + out[7][0][1])
331  - 0.5*(x - y - 2*z + 2)*(x - 1)*y;
332 
333  out[8][0][0] = -(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1)
334  + (x - z + 1)*((y + z + 1)*y + 2*z));
335  out[8][0][1] = -(x - z + 1)*(2*y + z + 1)*(x - 1);
336  out[8][0][2] = 0.5*(out[8][0][0] + out[8][0][1])
337  - 0.5*(-(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1))
338  + (x - z + 1)*(((x - 1)*y - 1) + 2*x + 1));
339 
340  // upper edges
341  out[9][0][0] = 2*z*(y + z - 1);
342  out[9][0][1] = 2*z*(x - z - 1);
343  out[9][0][2] = 0.5*(out[9][0][0] + out[9][0][1])
344  + (y + z - 1)*(x - z - 1) + z*(x - y - 2*z);
345 
346  out[10][0][0] = -2*z*(y + z - 1);
347  out[10][0][1] = -2*z*(x - z + 1);
348  out[10][0][2] = 0.5*(out[10][0][0] + out[10][0][1])
349  - (x - z + 1)*(y + z - 1) - z*(x - y - 2*z + 2);
350 
351  out[11][0][0] = -2*z*(y + z + 1);
352  out[11][0][1] = -2*z*(x - z - 1);
353  out[11][0][2] = 0.5*(out[11][0][0] + out[11][0][1])
354  - ((y + z + 1)*(x - z - 1) + 4*z) - z*(x - y - 2*z + 2);
355 
356  out[12][0][0] = 2*z*(y + z + 1);
357  out[12][0][1] = 2*z*(x - z + 1);
358  out[12][0][2] = 0.5*(out[12][0][0] + out[12][0][1])
359  + (x - z + 1)*(y + z + 1) + z*(x - y - 2*z);
360 
361  // base face
362  out[13][0][0] = 2*((y + z - 1)*((y + 1)*(x - 1) - z*(x - y - z - 1))
363  + (x - z + 1)*(y + z - 1)*(y + 1 - z));
364  out[13][0][1] = 2*((x - z + 1)*((y + 1)*(x - 1) - z*(x - y - z - 1))
365  + (x - z + 1)*(y + z - 1)*(x - 1 + z));
366  out[13][0][2] = 0.5*(out[13][0][0] + out[13][0][1])
367  + (x - y - 2*z + 2)*((y + 1)*(x - 1) - z*(x - y - z - 1))
368  + (x - z + 1)*(y + z - 1)*(-(x - y - 2*z - 1));
369  }
370 
371  return;
372  }
373 
374  DUNE_THROW(NotImplemented, "LagrangePyramidLocalBasis::evaluateJacobian for order " << k);
375  }
376 
383  void partial(const std::array<unsigned int,3>& order,
384  const typename Traits::DomainType& in,
385  std::vector<typename Traits::RangeType>& out) const
386  {
387  auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
388 
389  out.resize(size());
390 
391  if (totalOrder == 0)
392  {
393  evaluateFunction(in, out);
394  return;
395  }
396 
397  if (k==0)
398  {
399  out[0] = 0;
400  return;
401  }
402 
403  if (k==1)
404  {
405  if (totalOrder == 1)
406  {
407  auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
408  if (in[0] > in[1])
409  {
410  switch (direction)
411  {
412  case 0:
413  out[0] = -1 + in[1];
414  out[1] = 1 - in[1];
415  out[2] = -in[1];
416  out[3] = in[1];
417  out[4] = 0;
418  break;
419  case 1:
420  out[0] = -1 + in[0] + in[2];
421  out[1] = -in[0] - in[2];
422  out[2] = 1 - in[0] - in[2];
423  out[3] = in[0]+in[2];
424  out[4] = 0;
425  break;
426  case 2:
427  out[0] = -1 + in[1];
428  out[1] = -in[1];
429  out[2] = -in[1];
430  out[3] = in[1];
431  out[4] = 1;
432  break;
433  default:
434  DUNE_THROW(RangeError, "Component out of range.");
435  }
436  }
437  else /* (in[0] <= in[1]) */
438  {
439  switch (direction)
440  {
441  case 0:
442  out[0] = -1 + in[1] + in[2];
443  out[1] = 1 - in[1] - in[2];
444  out[2] = -in[1] - in[2];
445  out[3] = in[1] + in[2];
446  out[4] = 0;
447  break;
448  case 1:
449  out[0] = -1 + in[0];
450  out[1] = -in[0];
451  out[2] = 1 - in[0];
452  out[3] = in[0];
453  out[4] = 0;
454  break;
455  case 2:
456  out[0] = -1 + in[0];
457  out[1] = -in[0];
458  out[2] = -in[0];
459  out[3] = in[0];
460  out[4] = 1;
461  break;
462  default:
463  DUNE_THROW(RangeError, "Component out of range.");
464  }
465  }
466  } else if (totalOrder == 2)
467  {
468  if ((order[0] == 1 && order[1] == 1) ||
469  (order[1] == 1 && order[2] == 1 && in[0] > in[1]) ||
470  (order[0] == 1 && order[2] == 1 && in[0] <=in[1]))
471  {
472  out = {1, -1, -1, 1, 0};
473  } else
474  {
475  out = {0, 0, 0, 0, 0};
476  }
477 
478  } else
479  {
480  out = {0, 0, 0, 0, 0};
481  }
482 
483  return;
484  }
485 
486  if (k==2)
487  {
488  if (totalOrder == 1)
489  {
490  // transform to reference element with base [-1,1]^2
491  const R x = 2.0*in[0] + in[2] - 1.0;
492  const R y = 2.0*in[1] + in[2] - 1.0;
493  const R z = in[2];
494 
495  auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
496 
497  // transformation of the gradient leads to a multiplication
498  // with the Jacobian [2 0 0; 0 2 0; 1 1 1]
499  if (x > y)
500  {
501  switch (direction)
502  {
503  case 0:
504  out[0] = 0.5*(y - z - 1)*(y - z)*(2*x + 2*z - 1);
505  out[1] = 2*(-0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z) + (x + z)*(y - z)*(-y + z + 1)) - z);
506  out[2] = 0.5*(y - z)*(y - z + 1)*(2*x + 2*z - 1);
507  out[3] = 0.5*(y - z)*(2*x + 2*z + 1)*(y - z + 1);
508  out[4] = 0;
509  out[5] = -(y - z + 1)*(2*x + z - 1)*(y - 1);
510  out[6] = -(y - z + 1)*(2*x + z + 1)*(y - 1);
511  out[7] = -(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1) + (x + z - 1)*((y - z - 1)*y + 2*z));
512  out[8] = -(y - z + 1)*(2*x + z)*y;
513  out[9] = 2*z*(y - z - 1);
514  out[10] = -2*z*(y - z - 1);
515  out[11] = -2*z*(y - z + 1);
516  out[12] = 2*z*(y - z + 1);
517  out[13] = 2*((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*(y - 1 + z));
518  break;
519  case 1:
520  out[0] = 0.5*(x + z)*(x + z - 1)*(2*y - 2*z - 1);
521  out[1] = 2*(-0.25*((x + z)*((x + z + 1)*(-y + z + 1) - 4*z) + (x + z)*(y - z)*(-(x + z + 1))) + z);
522  out[2] = 0.5*(x + z)*(2*y - 2*z + 1)*(x + z - 1);
523  out[3] = 0.5*(2*y - 2*z + 1)*(x + z)*(x + z + 1);
524  out[4] = 0;
525  out[5] = -(x + z - 1)*(y - 1)*x - (y - z + 1)*(x + z - 1)*x;
526  out[6] = -(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1) + (y - z + 1)*((x + z + 1)*x + 2*z));
527  out[7] = -(x + z - 1)*(2*y - z - 1)*(x + 1);
528  out[8] = -(2*y - z + 1)*(x + z - 1)*(x + 1);
529  out[9] = 2*z*(x + z - 1);
530  out[10] = -2*z*(x + z + 1);
531  out[11] = -2*z*(x + z - 1);
532  out[12] = 2*z*(x + z + 1);
533  out[13] = 2*((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*(x + 1 - z));
534  break;
535  case 2:
536  out[0] = -((y - z)*(2*x + 2*z - 1)*(z - y + 1))/2;
537  out[1] = ((y - z + 1)*(y - 2*x + z + 2*x*y - 2*x*z + 2*y*z - 2*z*z))/2;
538  out[2] = ((y - z)*(2*x + 2*z - 1)*(y - z + 1))/2;
539  out[3] = ((y - z)*(2*x + 2*z + 1)*(y - z + 1))/2;
540  out[4] = 4*z - 1;
541  out[5] = (-(y - z + 1)*(2*x + z - 1)*(y - 1) - (x + z - 1)*(y - 1)*x - (y - z + 1)*(x + z - 1)*x + (x + z - 1)*(y - 1)*x - (y - z + 1)*(y - 1)*x)/2;
542  out[6] = -((y - z + 1)*(3*y - 2*x + z + 3*x*y + x*z + y*z + x*x - 1))/2;
543  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);
544  out[8] = -((y - z + 1)*(y + z + 3*x*y + x*z + y*z + x*x - 1))/2;
545  out[9] = -(x + 3*z - 1)*(z - y + 1);
546  out[10] = (x + z + 1)*(z - y + 1) - 2*y*z - 6*z + 2*z*z;
547  out[11] = -(x + 3*z - 1)*(y - z + 1);
548  out[12] = (x + 3*z + 1)*(y - z + 1);
549  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);
550  break;
551  default:
552  DUNE_THROW(RangeError, "Component out of range.");
553  }
554  }
555  else // x <= y
556  {
557  switch (direction)
558  {
559  case 0:
560  out[0] = -((y + z)*(2*z - 2*x + 1)*(y + z - 1))/2;
561  out[1] = ((y + z)*(2*x - 2*z + 1)*(y + z - 1))/2;
562  out[2] = -((y + z + 1)*(y - 3*z - 2*x*y - 2*x*z + 2*y*z + 2*z*z))/2;
563  out[3] = ((y + z)*(2*x - 2*z + 1)*(y + z + 1))/2;
564  out[4] = 0;
565  out[5] = (y + 1)*(y + z - 1)*(z - 2*x + 1);
566  out[6] = -(y + 1)*(2*x - z + 1)*(y + z - 1);
567  out[7] = -y*(2*x - z)*(y + z - 1);
568  out[8] = z - z*(2*x + 1) - (2*z + y*(y + z + 1))*(x - z + 1) - y*(x - 1)*(y + z + 1);
569  out[9] = 2*z*(y + z - 1);
570  out[10] = -2*z*(y + z - 1);
571  out[11] = -2*z*(y + z + 1);
572  out[12] = 2*z*(y + z + 1);
573  out[13] = 2*(y + z - 1)*(2*x - z + 2*x*y - 2*x*z + 2*z*z);
574  break;
575  case 1:
576  out[0] = -(x - z)*(y + z - 0.5)*(z - x + 1);
577  out[1] = ((x - z)*(2*y + 2*z - 1)*(x - z + 1))/2;
578  out[2] = -((z - x + 1)*(x + 3*z + 2*x*y + 2*x*z - 2*y*z - 2*z*z))/2;
579  out[3] = ((x - z)*(2*y + 2*z + 1)*(x - z + 1))/2;
580  out[4] = 0;
581  out[5] = z - z*(2*y + 1) - (2*z - x*(z - x + 1))*(y + z - 1) + x*(y + 1)*(z - x + 1);
582  out[6] = -x*(2*y + z)*(x - z + 1);
583  out[7] = -(x - 1)*(x - z + 1)*(2*y + z - 1);
584  out[8] = -(x - 1)*(x - z + 1)*(2*y + z + 1);
585  out[9] = -2*z*(z - x + 1);
586  out[10] = -2*z*(x - z + 1);
587  out[11] = 2*z*(z - x + 1);
588  out[12] = 2*z*(x - z + 1);
589  out[13] = 2*(x - z + 1)*(2*x*y - z - 2*y + 2*y*z + 2*z*z);
590  break;
591  case 2:
592  out[0] = -((x - z)*(2*y + 2*z - 1)*(z - x + 1))/2;
593  out[1] = ((x - z)*(2*y + 2*z - 1)*(x - z + 1))/2;
594  out[2] = ((x - z + 1)*(x - 2*y + z + 2*x*y + 2*x*z - 2*y*z - 2*z*z))/2;
595  out[3] = ((x - z)*(2*y + 2*z + 1)*(x - z + 1))/2;
596  out[4] = 4*z - 1;
597  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);
598  out[6] = -((x - z + 1)*(x + z + 3*x*y + x*z + y*z + y*y - 1))/2;
599  out[7] = -((x - z + 1)*(3*x*y - 4*y - z - x + x*z + y*z + y*y + 1))/2;
600  out[8] = -((x - z + 1)*(3*x - 2*y + z + 3*x*y + x*z + y*z + y*y - 1))/2;
601  out[9] = -(z - x + 1)*(y + 3*z - 1);
602  out[10] = -(x - z + 1)*(y + 3*z - 1);
603  out[11] = (y + z + 1)*(z - x + 1) - 2*x*z - 6*z + 2*z*z;
604  out[12] = (x - z + 1)*(y + 3*z + 1);
605  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);
606  break;
607  default:
608  DUNE_THROW(RangeError, "Component out of range.");
609  }
610  }
611  } else {
612  DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
613  }
614 
615  return;
616  }
617 
618  DUNE_THROW(NotImplemented, "LagrangePyramidLocalBasis::partial for order " << k);
619  }
620 
622  static constexpr unsigned int order ()
623  {
624  return k;
625  }
626  };
627 
632  template<unsigned int k>
633  class LagrangePyramidLocalCoefficients
634  {
635  public:
637  LagrangePyramidLocalCoefficients ()
638  : localKeys_(size())
639  {
640  if (k==0)
641  {
642  localKeys_[0] = LocalKey(0,0,0);
643  return;
644  }
645 
646  if (k==1)
647  {
648  for (std::size_t i=0; i<size(); i++)
649  localKeys_[i] = LocalKey(i,3,0);
650  return;
651  }
652 
653  if (k==2)
654  {
655  // Vertex shape functions
656  localKeys_[0] = LocalKey(0,3,0);
657  localKeys_[1] = LocalKey(1,3,0);
658  localKeys_[2] = LocalKey(2,3,0);
659  localKeys_[3] = LocalKey(3,3,0);
660  localKeys_[4] = LocalKey(4,3,0);
661 
662  // Edge shape functions
663  localKeys_[5] = LocalKey(0,2,0);
664  localKeys_[6] = LocalKey(1,2,0);
665  localKeys_[7] = LocalKey(2,2,0);
666  localKeys_[8] = LocalKey(3,2,0);
667  localKeys_[9] = LocalKey(4,2,0);
668  localKeys_[10] = LocalKey(5,2,0);
669  localKeys_[11] = LocalKey(6,2,0);
670  localKeys_[12] = LocalKey(7,2,0);
671 
672  // base face shape function
673  localKeys_[13] = LocalKey(0,1,0);
674 
675  return;
676  }
677 
678  // No general case
679  DUNE_THROW(NotImplemented, "LagrangePyramidLocalCoefficients for order " << k);
680 
681  }
682 
684  static constexpr std::size_t size ()
685  {
686  std::size_t result = 0;
687  for (unsigned int i=0; i<=k; i++)
688  result += power(i+1,2);
689  return result;
690  }
691 
693  const LocalKey& localKey (std::size_t i) const
694  {
695  return localKeys_[i];
696  }
697 
698  private:
699  std::vector<LocalKey> localKeys_;
700  };
701 
706  template<class LocalBasis>
707  class LagrangePyramidLocalInterpolation
708  {
709  public:
710 
718  template<typename F, typename C>
719  void interpolate (const F& ff, std::vector<C>& out) const
720  {
721  constexpr auto k = LocalBasis::order();
722  using D = typename LocalBasis::Traits::DomainType;
723  using DF = typename LocalBasis::Traits::DomainFieldType;
724 
725  auto&& f = Impl::makeFunctionWithCallOperator<D>(ff);
726 
727  out.resize(LocalBasis::size());
728 
729  // Specialization for zero-order case
730  if (k==0)
731  {
732  auto center = ReferenceElements<DF,3>::general(GeometryTypes::pyramid).position(0,0);
733  out[0] = f(center);
734  return;
735  }
736 
737  // Specialization for first-order case
738  if (k==1)
739  {
740  for (unsigned int i=0; i<LocalBasis::size(); i++)
741  {
742  auto vertex = ReferenceElements<DF,3>::general(GeometryTypes::pyramid).position(i,3);
743  out[i] = f(vertex);
744  }
745  return;
746  }
747 
748  // Specialization for second-order case
749  if (k==2)
750  {
751  out[0] = f( D( {0.0, 0.0, 0.0} ) );
752  out[1] = f( D( {1.0, 0.0, 0.0} ) );
753  out[2] = f( D( {0.0, 1.0, 0.0} ) );
754  out[3] = f( D( {1.0, 1.0, 0.0} ) );
755  out[4] = f( D( {0.0, 0.0, 1.0} ) );
756  out[5] = f( D( {0.0, 0.5, 0.0} ) );
757  out[6] = f( D( {1.0, 0.5, 0.0} ) );
758  out[7] = f( D( {0.5, 0.0, 0.0} ) );
759  out[8] = f( D( {0.5, 1.0, 0.0} ) );
760  out[9] = f( D( {0.0, 0.0, 0.5} ) );
761  out[10] = f( D( {0.5, 0.0, 0.5} ) );
762  out[11] = f( D( {0.0, 0.5, 0.5} ) );
763  out[12] = f( D( {0.5, 0.5, 0.5} ) );
764  out[13] = f( D( {0.5, 0.5, 0.0} ) );
765 
766  return;
767  }
768 
769  DUNE_THROW(NotImplemented, "LagrangePyramidLocalInterpolation not implemented for order " << k);
770  }
771 
772  };
773 
774 } } // namespace Dune::Impl
775 
776 namespace Dune
777 {
808  template<class D, class R, int k>
809  class LagrangePyramidLocalFiniteElement
810  {
811  public:
814  using Traits = LocalFiniteElementTraits<Impl::LagrangePyramidLocalBasis<D,R,k>,
815  Impl::LagrangePyramidLocalCoefficients<k>,
816  Impl::LagrangePyramidLocalInterpolation<Impl::LagrangePyramidLocalBasis<D,R,k> > >;
817 
823  LagrangePyramidLocalFiniteElement() {}
824 
827  const typename Traits::LocalBasisType& localBasis () const
828  {
829  return basis_;
830  }
831 
834  const typename Traits::LocalCoefficientsType& localCoefficients () const
835  {
836  return coefficients_;
837  }
838 
841  const typename Traits::LocalInterpolationType& localInterpolation () const
842  {
843  return interpolation_;
844  }
845 
847  static constexpr std::size_t size ()
848  {
849  return Impl::LagrangePyramidLocalBasis<D,R,k>::size();
850  }
851 
854  static constexpr GeometryType type ()
855  {
856  return GeometryTypes::pyramid;
857  }
858 
859  private:
860  Impl::LagrangePyramidLocalBasis<D,R,k> basis_;
861  Impl::LagrangePyramidLocalCoefficients<k> coefficients_;
862  Impl::LagrangePyramidLocalInterpolation<Impl::LagrangePyramidLocalBasis<D,R,k> > interpolation_;
863  };
864 
865 } // namespace Dune
866 
867 #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 compile-time polynomial order.
Definition: lagrangepyramid.hh:810
Dune::LagrangePyramidLocalFiniteElement::localCoefficients
const Traits::LocalCoefficientsType & localCoefficients() const
Returns the assignment of the degrees of freedom to the element subentities.
Definition: lagrangepyramid.hh:834
Dune::LagrangePyramidLocalFiniteElement::size
static constexpr std::size_t size()
The number of shape functions.
Definition: lagrangepyramid.hh:847
Dune::LagrangePyramidLocalFiniteElement::LagrangePyramidLocalFiniteElement
LagrangePyramidLocalFiniteElement()
Default constructor.
Definition: lagrangepyramid.hh:823
Dune::LagrangePyramidLocalFiniteElement::type
static constexpr GeometryType type()
The reference element that the local finite element is defined on.
Definition: lagrangepyramid.hh:854
Dune::LagrangePyramidLocalFiniteElement::localBasis
const Traits::LocalBasisType & localBasis() const
Returns the local basis, i.e., the set of shape functions.
Definition: lagrangepyramid.hh:827
Dune::LagrangePyramidLocalFiniteElement::localInterpolation
const Traits::LocalInterpolationType & localInterpolation() const
Returns object that evaluates degrees of freedom.
Definition: lagrangepyramid.hh:841
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)