Dune Core Modules (2.6.0)

pyramidp2localbasis.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_PYRAMID_P2_LOCALBASIS_HH
4 #define DUNE_PYRAMID_P2_LOCALBASIS_HH
5 
6 #include <numeric>
7 
8 #include <dune/common/fmatrix.hh>
9 
10 #include <dune/localfunctions/common/localbasis.hh>
11 
12 namespace Dune
13 {
28  template<class D, class R>
29  class PyramidP2LocalBasis
30  {
31  public:
33  typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,R,1,Dune::FieldVector<R,1>,
34  Dune::FieldMatrix<R,1,3> > Traits;
35 
37  unsigned int size () const
38  {
39  return 14;
40  }
41 
43  inline void evaluateFunction (const typename Traits::DomainType& in,
44  std::vector<typename Traits::RangeType>& out) const
45  {
46  out.resize(14);
47 
48  // transform to reference element with base [-1,1]^2
49  const R x = 2.0*in[0] + in[2] - 1.0;
50  const R y = 2.0*in[1] + in[2] - 1.0;
51  const R z = in[2];
52 
53  if (x > y)
54  {
55  // vertices
56  out[0] = 0.25*(x + z)*(x + z - 1)*(y - z - 1)*(y - z);
57  out[1] = -0.25*(x + z)*(y - z)*((x + z + 1)*(-y + z + 1) - 4*z) - z*(x - y);
58  out[2] = 0.25*(x + z)*(y - z)*(y - z + 1)*(x + z - 1);
59  out[3] = 0.25*(y - z)*(x + z)*(y - z + 1)*(x + z + 1);
60  out[4] = z*(2*z - 1);
61 
62  // lower edges
63  out[5] = -0.5*(y - z + 1)*(x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1));
64  out[6] = -0.5*(y - z + 1)*(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1));
65  out[7] = -0.5*(x + z - 1)*(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1));
66  out[8] = -0.5*(y - z + 1)*(x + z - 1)*(x + 1)*y;
67 
68  // upper edges
69  out[9] = z*(x + z - 1)*(y - z - 1);
70  out[10] = -z*((x + z + 1)*(y - z - 1) + 4*z);
71  out[11] = -z*(y - z + 1)*(x + z - 1);
72  out[12] = z*(y - z + 1)*(x + z + 1);
73 
74  // base face
75  out[13] = (y - z + 1)*(x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1));
76  }
77  else
78  {
79  // vertices
80  out[0] = 0.25*(y + z)*(y + z - 1)*(x - z - 1)*(x - z);
81  out[1] = -0.25*(x - z)*(y + z)*(x - z + 1)*(-y - z + 1);
82  out[2] = 0.25*(x - z)*(y + z)*((x - z - 1)*(y + z + 1) + 4*z) + z*(x - y);
83  out[3] = 0.25*(y + z)*(x - z)*(x - z + 1)*(y + z + 1);
84  out[4] = z*(2*z - 1);
85 
86  // lower edges
87  out[5] = -0.5*(y + z - 1)*(((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1));
88  out[6] = -0.5*(x - z + 1)*(y + z - 1)*(y + 1)*x;
89  out[7] = -0.5*(x - z + 1)*(y + z - 1)*(x - 1)*y;
90  out[8] = -0.5*(x - z + 1)*(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1));
91 
92  // upper edges
93  out[9] = z*(y + z - 1)*(x - z - 1);
94  out[10] = -z*(x - z + 1)*(y + z - 1);
95  out[11] = -z*((y + z + 1)*(x - z - 1) + 4*z);
96  out[12] = z*(x - z + 1)*(y + z + 1);
97 
98  // base face
99  out[13] = (x - z + 1)*(y + z - 1)*((y + 1)*(x - 1) - z*(x - y - z - 1));
100  }
101  }
102 
104  inline void
105  evaluateJacobian (const typename Traits::DomainType& in, // position
106  std::vector<typename Traits::JacobianType>& out) const // return value
107  {
108  out.resize(14);
109 
110  // transform to reference element with base [-1,1]^2
111  const R x = 2.0*in[0] + in[2] - 1.0;
112  const R y = 2.0*in[1] + in[2] - 1.0;
113  const R z = in[2];
114 
115  // transformation of the gradient leads to a multiplication
116  // with the Jacobian [2 0 0; 0 2 0; 1 1 1]
117  if (x > y)
118  {
119  // vertices
120  out[0][0][0] = 0.5*(y - z - 1)*(y - z)*(2*x + 2*z - 1);
121  out[0][0][1] = 0.5*(x + z)*(x + z - 1)*(2*y - 2*z - 1);
122  out[0][0][2] = 0.5*(out[0][0][0] + out[0][0][1])
123  + 0.25*((2*x + 2*z - 1)*(y - z - 1)*(y - z)
124  + (x + z)*(x + z - 1)*(-2*y + 2*z + 1));
125 
126  out[1][0][0] = 2*(-0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z)
127  + (x + z)*(y - z)*(-y + z + 1)) - z);
128  out[1][0][1] = 2*(-0.25*((x + z)*((x + z + 1)*(-y + z + 1) - 4*z)
129  + (x + z)*(y - z)*(-(x + z + 1))) + z);
130  out[1][0][2] = 0.5*(out[1][0][0] + out[1][0][1])
131  - 0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z)
132  - (x + z)*((x + z + 1)*(-y + z + 1) - 4*z)
133  + (x + z)*(y - z)*(x - y + 2*z - 2))
134  - (x - y);
135 
136  out[2][0][0] = 0.5*(y - z)*(y - z + 1)*(2*x + 2*z - 1);
137  out[2][0][1] = 0.5*(x + z)*(2*y - 2*z + 1)*(x + z - 1);
138  out[2][0][2] = 0.5*(out[2][0][0] + out[2][0][1])
139  + 0.25*((y - x - 2*z)*(y - z + 1)*(x + z - 1)
140  + (x + z)*(y - z)*(y - x - 2*z + 2));
141 
142  out[3][0][0] = 0.5*(y - z)*(2*x + 2*z + 1)*(y - z + 1);
143  out[3][0][1] = 0.5*(2*y - 2*z + 1)*(x + z)*(x + z + 1);
144  out[3][0][2] = 0.5*(out[3][0][0] + out[3][0][1])
145  + 0.25*((y - x - 2*z)*(y - z + 1)*(x + z + 1)
146  + (y - z)*(x + z)*(y - x - 2*z));
147 
148  out[4][0][0] = 0;
149  out[4][0][1] = 0;
150  out[4][0][2] = 4*z - 1;
151 
152  // lower edges
153  out[5][0][0] = -((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
154  + (y - z + 1)*(x + z - 1)*((y - 1) + z));
155  out[5][0][1] = -((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
156  + (y - z + 1)*(x + z - 1)*((x + 1) - z));
157  out[5][0][2] = 0.5*(out[5][0][0] + out[5][0][1])
158  - 0.5*((-x + y - 2*z + 2)*((y - 1)*(x + 1) + z*(x - y + z + 1))
159  + (y - z + 1)*(x + z - 1)*(x - y + 2*z + 1));
160 
161  out[6][0][0] = -(y - z + 1)*(2*x + z + 1)*(y - 1);
162  out[6][0][1] = -(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1)
163  + (y - z + 1)*((x + z + 1)*x + 2*z));
164  out[6][0][2] = 0.5*(out[6][0][0] + out[6][0][1])
165  - 0.5*(-(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1))
166  + (y - z + 1)*(((y - 1)*x - 1) + 2*y + 1));
167 
168  out[7][0][0] = -(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1)
169  + (x + z - 1)*((y - z - 1)*y + 2*z));
170  out[7][0][1] = -(x + z - 1)*(2*y - z - 1)*(x + 1);
171  out[7][0][2] = 0.5*(out[7][0][0] + out[7][0][1])
172  - 0.5*(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1)
173  + (x + z - 1)*((-(x + 1)*y - 1) + 2*x + 1));
174 
175  out[8][0][0] = -(y - z + 1)*(2*x + z)*y;
176  out[8][0][1] = -(2*y - z + 1)*(x + z - 1)*(x + 1);
177  out[8][0][2] = 0.5*(out[8][0][0] + out[8][0][1])
178  - 0.5*(-x + y - 2*z + 2)*(x + 1)*y;
179 
180  // upper edges
181  out[9][0][0] = 2*z*(y - z - 1);
182  out[9][0][1] = 2*z*(x + z - 1);
183  out[9][0][2] = 0.5*(out[9][0][0] + out[9][0][1])
184  + (x + z - 1)*(y - z - 1) + z*(-x + y - 2*z);
185 
186  out[10][0][0] = -2*z*(y - z - 1);
187  out[10][0][1] = -2*z*(x + z + 1);
188  out[10][0][2] = 0.5*(out[10][0][0] + out[10][0][1])
189  - ((x + z + 1)*(y - z - 1) + 4*z)
190  - z*(-x + y - 2*z + 2);
191 
192  out[11][0][0] = -2*z*(y - z + 1);
193  out[11][0][1] = -2*z*(x + z - 1);
194  out[11][0][2] = 0.5*(out[11][0][0] + out[11][0][1])
195  - (y - z + 1)*(x + z - 1) - z*(-x + y - 2*z + 2);
196 
197  out[12][0][0] = 2*z*(y - z + 1);
198  out[12][0][1] = 2*z*(x + z + 1);
199  out[12][0][2] = 0.5*(out[12][0][0] + out[12][0][1])
200  + (y - z + 1)*(x + z + 1) + z*(-x + y - 2*z);
201 
202  // base face
203  out[13][0][0] = 2*((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
204  + (y - z + 1)*(x + z - 1)*(y - 1 + z));
205  out[13][0][1] = 2*((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
206  + (y - z + 1)*(x + z - 1)*(x + 1 - z));
207  out[13][0][2] = 0.5*(out[13][0][0] + out[13][0][1])
208  + ((-x + y - 2*z + 2)*((y - 1)*(x + 1) + z*(x - y + z + 1))
209  + (y - z + 1)*(x + z - 1)*(x - y + 2*z + 1));
210  }
211  else
212  {
213  // vertices
214  out[0][0][0] = 0.5*(y + z)*(y + z - 1)*(2*x - 2*z - 1);
215  out[0][0][1] = 0.5*(2*y + 2*z - 1)*(x - z - 1)*(x - z);
216  out[0][0][2] = 0.5*(out[0][0][0] + out[0][0][1])
217  + 0.25*((2*y + 2*z - 1)*(x - z - 1)*(x - z)
218  + (y + z)*(y + z - 1)*(-2*x + 2*z + 1));
219 
220  out[1][0][0] = -0.5*(y + z)*(2*x - 2*z + 1)*(-y - z + 1);
221  out[1][0][1] = -0.5*(x - z)*(x - z + 1)*(-2*y - 2*z + 1);
222  out[1][0][2] = 0.5*(out[1][0][0] + out[1][0][1])
223  - 0.25*((x - y - 2*z)*(x - z + 1)*(-y - z + 1)
224  + (x - z)*(y + z)*(-x + y + 2*z - 2));
225 
226  out[2][0][0] = 0.5*((y + z)*((x - z - 1)*(y + z + 1) + 4*z)
227  + (x - z)*(y + z)*(y + z + 1) + 4*z);
228  out[2][0][1] = 0.5*((x - z)*((x - z - 1)*(y + z + 1) + 4*z)
229  + (x - z)*(y + z)*(x - z - 1) - 4*z);
230  out[2][0][2] = 0.5*(out[2][0][0] + out[2][0][1])
231  + 0.25*((x - y - 2*z)*((x - z - 1)*(y + z + 1) + 4*z)
232  + (x - z)*(y + z)*(x - y - 2*z + 2) + 4*(x - y));
233 
234  out[3][0][0] = 0.5*(y + z)*(2*x - 2*z + 1)*(y + z + 1);
235  out[3][0][1] = 0.5*(x - z)*(x - z + 1)*(2*y + 2*z + 1);
236  out[3][0][2] = 0.5*(out[3][0][0] + out[3][0][1])
237  + 0.25*((x - y - 2*z)*(x - z + 1)*(y + z + 1)
238  + (y + z)*(x - z)*(x - y - 2*z));
239 
240  out[4][0][0] = 0;
241  out[4][0][1] = 0;
242  out[4][0][2] = 4*z - 1;
243 
244  // lower edges
245  out[5][0][0] = -(y + z - 1)*(2*x - z - 1)*(y + 1);
246  out[5][0][1] = -(((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1)
247  + (y + z - 1)*((x - z - 1)*x + 2*z));
248  out[5][0][2] = 0.5*(out[5][0][0] + out[5][0][1])
249  - 0.5*((((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1))
250  + (y + z - 1)*((-(y + 1)*x - 1) + 2*y + 1));
251 
252  out[6][0][0] = -(2*x - z + 1)*(y + z - 1)*(y + 1);
253  out[6][0][1] = -(x - z + 1)*(2*y + z)*x;
254  out[6][0][2] = 0.5*(out[6][0][0] + out[6][0][1])
255  - 0.5*(x - y - 2*z + 2)*(y + 1)*x;
256 
257  out[7][0][0] = -(2*x - z)*(y + z - 1)*y;
258  out[7][0][1] = -(x - z + 1)*(2*y + z - 1)*(x - 1);
259  out[7][0][2] = 0.5*(out[7][0][0] + out[7][0][1])
260  - 0.5*(x - y - 2*z + 2)*(x - 1)*y;
261 
262  out[8][0][0] = -(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1)
263  + (x - z + 1)*((y + z + 1)*y + 2*z));
264  out[8][0][1] = -(x - z + 1)*(2*y + z + 1)*(x - 1);
265  out[8][0][2] = 0.5*(out[8][0][0] + out[8][0][1])
266  - 0.5*(-(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1))
267  + (x - z + 1)*(((x - 1)*y - 1) + 2*x + 1));
268 
269  // upper edges
270  out[9][0][0] = 2*z*(y + z - 1);
271  out[9][0][1] = 2*z*(x - z - 1);
272  out[9][0][2] = 0.5*(out[9][0][0] + out[9][0][1])
273  + (y + z - 1)*(x - z - 1) + z*(x - y - 2*z);
274 
275  out[10][0][0] = -2*z*(y + z - 1);
276  out[10][0][1] = -2*z*(x - z + 1);
277  out[10][0][2] = 0.5*(out[10][0][0] + out[10][0][1])
278  - (x - z + 1)*(y + z - 1) - z*(x - y - 2*z + 2);
279 
280  out[11][0][0] = -2*z*(y + z + 1);
281  out[11][0][1] = -2*z*(x - z - 1);
282  out[11][0][2] = 0.5*(out[11][0][0] + out[11][0][1])
283  - ((y + z + 1)*(x - z - 1) + 4*z) - z*(x - y - 2*z + 2);
284 
285  out[12][0][0] = 2*z*(y + z + 1);
286  out[12][0][1] = 2*z*(x - z + 1);
287  out[12][0][2] = 0.5*(out[12][0][0] + out[12][0][1])
288  + (x - z + 1)*(y + z + 1) + z*(x - y - 2*z);
289 
290  // base face
291  out[13][0][0] = 2*((y + z - 1)*((y + 1)*(x - 1) - z*(x - y - z - 1))
292  + (x - z + 1)*(y + z - 1)*(y + 1 - z));
293  out[13][0][1] = 2*((x - z + 1)*((y + 1)*(x - 1) - z*(x - y - z - 1))
294  + (x - z + 1)*(y + z - 1)*(x - 1 + z));
295  out[13][0][2] = 0.5*(out[13][0][0] + out[13][0][1])
296  + (x - y - 2*z + 2)*((y + 1)*(x - 1) - z*(x - y - z - 1))
297  + (x - z + 1)*(y + z - 1)*(-(x - y - 2*z - 1));
298  }
299  }
300 
302  void partial (const std::array<unsigned int, 3>& order,
303  const typename Traits::DomainType& in, // position
304  std::vector<typename Traits::RangeType>& out) const // return value
305  {
306  auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
307  if (totalOrder == 0) {
308  evaluateFunction(in, out);
309  } else if (totalOrder == 1) {
310  out.resize(size());
311 
312  // transform to reference element with base [-1,1]^2
313  const R x = 2.0*in[0] + in[2] - 1.0;
314  const R y = 2.0*in[1] + in[2] - 1.0;
315  const R z = in[2];
316 
317  auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
318 
319  // transformation of the gradient leads to a multiplication
320  // with the Jacobian [2 0 0; 0 2 0; 1 1 1]
321  if (x > y) {
322  switch (direction) {
323  case 0:
324  out[0] = 0.5*(y - z - 1)*(y - z)*(2*x + 2*z - 1);
325  out[1] = 2*(-0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z) + (x + z)*(y - z)*(-y + z + 1)) - z);
326  out[2] = 0.5*(y - z)*(y - z + 1)*(2*x + 2*z - 1);
327  out[3] = 0.5*(y - z)*(2*x + 2*z + 1)*(y - z + 1);
328  out[4] = 0;
329  out[5] = -((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*((y - 1) + z));
330  out[6] = -(y - z + 1)*(2*x + z + 1)*(y - 1);
331  out[7] = -(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1) + (x + z - 1)*((y - z - 1)*y + 2*z));
332  out[8] = -(y - z + 1)*(2*x + z)*y;
333  out[9] = 2*z*(y - z - 1);
334  out[10] = -2*z*(y - z - 1);
335  out[11] = -2*z*(y - z + 1);
336  out[12] = 2*z*(y - z + 1);
337  out[13] = 2*((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*(y - 1 + z));
338  break;
339  case 1:
340  out[0] = 0.5*(x + z)*(x + z - 1)*(2*y - 2*z - 1);
341  out[1] = 2*(-0.25*((x + z)*((x + z + 1)*(-y + z + 1) - 4*z) + (x + z)*(y - z)*(-(x + z + 1))) + z);
342  out[2] = 0.5*(x + z)*(2*y - 2*z + 1)*(x + z - 1);
343  out[3] = 0.5*(2*y - 2*z + 1)*(x + z)*(x + z + 1);
344  out[4] = 0;
345  out[5] = -((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*((x + 1) - z));
346  out[6] = -(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1) + (y - z + 1)*((x + z + 1)*x + 2*z));
347  out[7] = -(x + z - 1)*(2*y - z - 1)*(x + 1);
348  out[8] = -(2*y - z + 1)*(x + z - 1)*(x + 1);
349  out[9] = 2*z*(x + z - 1);
350  out[10] = -2*z*(x + z + 1);
351  out[11] = -2*z*(x + z - 1);
352  out[12] = 2*z*(x + z + 1);
353  out[13] = 2*((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*(x + 1 - z));
354  break;
355  case 2:
356  out[0] = -((y - z)*(2*x + 2*z - 1)*(z - y + 1))/2;
357  out[1] = ((y - z + 1)*(y - 2*x + z + 2*x*y - 2*x*z + 2*y*z - 2*z*z))/2;
358  out[2] = ((y - z)*(2*x + 2*z - 1)*(y - z + 1))/2;
359  out[3] = ((y - z)*(2*x + 2*z + 1)*(y - z + 1))/2;
360  out[4] = 4*z - 1;
361  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;
362  out[6] = -((y - z + 1)*(3*y - 2*x + z + 3*x*y + x*z + y*z + x*x - 1))/2;
363  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);
364  out[8] = -((y - z + 1)*(y + z + 3*x*y + x*z + y*z + x*x - 1))/2;
365  out[9] = -(x + 3*z - 1)*(z - y + 1);
366  out[10] = (x + z + 1)*(z - y + 1) - 2*y*z - 6*z + 2*z*z;
367  out[11] = -(x + 3*z - 1)*(y - z + 1);
368  out[12] = (x + 3*z + 1)*(y - z + 1);
369  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);
370  break;
371  default:
372  DUNE_THROW(RangeError, "Component out of range.");
373  }
374  }
375  else // x <= y
376  {
377  switch (direction) {
378  case 0:
379  out[0] = -((y + z)*(2*z - 2*x + 1)*(y + z - 1))/2;
380  out[1] = ((y + z)*(2*x - 2*z + 1)*(y + z - 1))/2;
381  out[2] = -((y + z + 1)*(y - 3*z - 2*x*y - 2*x*z + 2*y*z + 2*z*z))/2;
382  out[3] = ((y + z)*(2*x - 2*z + 1)*(y + z + 1))/2;
383  out[4] = 0;
384  out[5] = (y + 1)*(y + z - 1)*(z - 2*x + 1);
385  out[6] = -(y + 1)*(2*x - z + 1)*(y + z - 1);
386  out[7] = -y*(2*x - z)*(y + z - 1);
387  out[8] = z - z*(2*x + 1) - (2*z + y*(y + z + 1))*(x - z + 1) - y*(x - 1)*(y + z + 1);
388  out[9] = 2*z*(y + z - 1);
389  out[10] = -2*z*(y + z - 1);
390  out[11] = -2*z*(y + z + 1);
391  out[12] = 2*z*(y + z + 1);
392  out[13] = 2*(y + z - 1)*(2*x - z + 2*x*y - 2*x*z + 2*z*z);
393  break;
394  case 1:
395  out[0] = -(x - z)*(y + z - 0.5)*(z - x + 1);
396  out[1] = ((x - z)*(2*y + 2*z - 1)*(x - z + 1))/2;
397  out[2] = -((z - x + 1)*(x + 3*z + 2*x*y + 2*x*z - 2*y*z - 2*z*z))/2;
398  out[3] = ((x - z)*(2*y + 2*z + 1)*(x - z + 1))/2;
399  out[4] = 0;
400  out[5] = z - z*(2*y + 1) - (2*z - x*(z - x + 1))*(y + z - 1) + x*(y + 1)*(z - x + 1);
401  out[6] = -x*(2*y + z)*(x - z + 1);
402  out[7] = -(x - 1)*(x - z + 1)*(2*y + z - 1);
403  out[8] = -(x - 1)*(x - z + 1)*(2*y + z + 1);
404  out[9] = -2*z*(z - x + 1);
405  out[10] = -2*z*(x - z + 1);
406  out[11] = 2*z*(z - x + 1);
407  out[12] = 2*z*(x - z + 1);
408  out[13] = 2*(x - z + 1)*(2*x*y - z - 2*y + 2*y*z + 2*z*z);
409  break;
410  case 2:
411  out[0] = -((x - z)*(2*y + 2*z - 1)*(z - x + 1))/2;
412  out[1] = ((x - z)*(2*y + 2*z - 1)*(x - z + 1))/2;
413  out[2] = ((x - z + 1)*(x - 2*y + z + 2*x*y + 2*x*z - 2*y*z - 2*z*z))/2;
414  out[3] = ((x - z)*(2*y + 2*z + 1)*(x - z + 1))/2;
415  out[4] = 4*z - 1;
416  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);
417  out[6] = -((x - z + 1)*(x + z + 3*x*y + x*z + y*z + y*y - 1))/2;
418  out[7] = -((x - z + 1)*(3*x*y - 4*y - z - x + x*z + y*z + y*y + 1))/2;
419  out[8] = -((x - z + 1)*(3*x - 2*y + z + 3*x*y + x*z + y*z + y*y - 1))/2;
420  out[9] = -(z - x + 1)*(y + 3*z - 1);
421  out[10] = -(x - z + 1)*(y + 3*z - 1);
422  out[11] = (y + z + 1)*(z - x + 1) - 2*x*z - 6*z + 2*z*z;
423  out[12] = (x - z + 1)*(y + 3*z + 1);
424  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);
425  break;
426  default:
427  DUNE_THROW(RangeError, "Component out of range.");
428  }
429  }
430  } else {
431  DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
432  }
433  }
434 
436  unsigned int order () const
437  {
438  return 2;
439  }
440  };
441 }
442 #endif
Dune::FieldMatrix
A dense n x m matrix.
Definition: fmatrix.hh:68
Dune::FieldVector
vector space out of a tensor product of fields.
Definition: fvector.hh:93
Dune::NotImplemented
Default exception for dummy implementations.
Definition: exceptions.hh:261
Dune::PyramidP2LocalBasis
Quadratic Lagrange shape functions on the pyramid.
Definition: pyramidp2localbasis.hh:30
Dune::PyramidP2LocalBasis::partial
void partial(const std::array< unsigned int, 3 > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivatives of all shape functions.
Definition: pyramidp2localbasis.hh:302
Dune::PyramidP2LocalBasis::order
unsigned int order() const
Polynomial order of the shape functions.
Definition: pyramidp2localbasis.hh:436
Dune::PyramidP2LocalBasis::size
unsigned int size() const
number of shape functions
Definition: pyramidp2localbasis.hh:37
Dune::PyramidP2LocalBasis::Traits
LocalBasisTraits< D, 3, Dune::FieldVector< D, 3 >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, 3 > > Traits
export type traits for function signature
Definition: pyramidp2localbasis.hh:34
Dune::PyramidP2LocalBasis::evaluateFunction
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: pyramidp2localbasis.hh:43
Dune::PyramidP2LocalBasis::evaluateJacobian
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: pyramidp2localbasis.hh:105
Dune::RangeError
Default exception class for range errors.
Definition: exceptions.hh:252
fmatrix.hh
Implements a matrix constructed from a given type representing a field and compile-time given number ...
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune::Hybrid::accumulate
T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:331
Dune
Dune namespace.
Definition: alignedallocator.hh:10
Dune::LocalBasisTraits
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:32
Dune::LocalBasisTraits::DomainType
D DomainType
domain type
Definition: localbasis.hh:43
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.80.0 (May 3, 22:32, 2024)