Dune Core Modules (unstable)

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