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