dune-localfunctions  2.2.1
pyramidp2localbasis.hh
Go to the documentation of this file.
1 // -*- tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=2 sw=2 sts=2:
3 #ifndef DUNE_PYRAMID_P2_LOCALBASIS_HH
4 #define DUNE_PYRAMID_P2_LOCALBASIS_HH
5 
6 #include <dune/common/fmatrix.hh>
7 
9 
10 namespace Dune
11 {
26  template<class D, class R>
28  {
29  public:
31  typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,R,1,Dune::FieldVector<R,1>,
32  Dune::FieldMatrix<R,1,3> > Traits;
33 
35  unsigned int size () const
36  {
37  return 14;
38  }
39 
41  inline void evaluateFunction (const typename Traits::DomainType& in,
42  std::vector<typename Traits::RangeType>& out) const
43  {
44  out.resize(14);
45 
46  // transform to reference element with base [-1,1]^2
47  const R x = 2.0*in[0] + in[2] - 1.0;
48  const R y = 2.0*in[1] + in[2] - 1.0;
49  const R z = in[2];
50 
51  if (x > y)
52  {
53  // vertices
54  out[0] = 0.25*(x + z)*(x + z - 1)*(y - z - 1)*(y - z);
55  out[1] = -0.25*(x + z)*(y - z)*((x + z + 1)*(-y + z + 1) - 4*z) - z*(x - y);
56  out[2] = 0.25*(x + z)*(y - z)*(y - z + 1)*(x + z - 1);
57  out[3] = 0.25*(y - z)*(x + z)*(y - z + 1)*(x + z + 1);
58  out[4] = z*(2*z - 1);
59 
60  // lower edges
61  out[5] = -0.5*(y - z + 1)*(x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1));
62  out[6] = -0.5*(y - z + 1)*(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1));
63  out[7] = -0.5*(x + z - 1)*(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1));
64  out[8] = -0.5*(y - z + 1)*(x + z - 1)*(x + 1)*y;
65 
66  // upper edges
67  out[9] = z*(x + z - 1)*(y - z - 1);
68  out[10] = -z*((x + z + 1)*(y - z - 1) + 4*z);
69  out[11] = -z*(y - z + 1)*(x + z - 1);
70  out[12] = z*(y - z + 1)*(x + z + 1);
71 
72  // base face
73  out[13] = (y - z + 1)*(x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1));
74  }
75  else
76  {
77  // vertices
78  out[0] = 0.25*(y + z)*(y + z - 1)*(x - z - 1)*(x - z);
79  out[1] = -0.25*(x - z)*(y + z)*(x - z + 1)*(-y - z + 1);
80  out[2] = 0.25*(x - z)*(y + z)*((x - z - 1)*(y + z + 1) + 4*z) + z*(x - y);
81  out[3] = 0.25*(y + z)*(x - z)*(x - z + 1)*(y + z + 1);
82  out[4] = z*(2*z - 1);
83 
84  // lower edges
85  out[5] = -0.5*(y + z - 1)*(((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1));
86  out[6] = -0.5*(x - z + 1)*(y + z - 1)*(y + 1)*x;
87  out[7] = -0.5*(x - z + 1)*(y + z - 1)*(x - 1)*y;
88  out[8] = -0.5*(x - z + 1)*(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1));
89 
90  // upper edges
91  out[9] = z*(y + z - 1)*(x - z - 1);
92  out[10] = -z*(x - z + 1)*(y + z - 1);
93  out[11] = -z*((y + z + 1)*(x - z - 1) + 4*z);
94  out[12] = z*(x - z + 1)*(y + z + 1);
95 
96  // base face
97  out[13] = (x - z + 1)*(y + z - 1)*((y + 1)*(x - 1) - z*(x - y - z - 1));
98  }
99  }
100 
102  inline void
103  evaluateJacobian (const typename Traits::DomainType& in, // position
104  std::vector<typename Traits::JacobianType>& out) const // return value
105  {
106  out.resize(14);
107 
108  // transform to reference element with base [-1,1]^2
109  const R x = 2.0*in[0] + in[2] - 1.0;
110  const R y = 2.0*in[1] + in[2] - 1.0;
111  const R z = in[2];
112 
113  // transformation of the gradient leads to a multiplication
114  // with the Jacobian [2 0 0; 0 2 0; 1 1 1]
115  if (x > y)
116  {
117  // vertices
118  out[0][0][0] = 0.5*(y - z - 1)*(y - z)*(2*x + 2*z - 1);
119  out[0][0][1] = 0.5*(x + z)*(x + z - 1)*(2*y - 2*z - 1);
120  out[0][0][2] = 0.5*(out[0][0][0] + out[0][0][1])
121  + 0.25*((2*x + 2*z - 1)*(y - z - 1)*(y - z)
122  + (x + z)*(x + z - 1)*(-2*y + 2*z + 1));
123 
124  out[1][0][0] = 2*(-0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z)
125  + (x + z)*(y - z)*(-y + z + 1)) - z);
126  out[1][0][1] = 2*(-0.25*((x + z)*((x + z + 1)*(-y + z + 1) - 4*z)
127  + (x + z)*(y - z)*(-(x + z + 1))) + z);
128  out[1][0][2] = 0.5*(out[1][0][0] + out[1][0][1])
129  - 0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z)
130  - (x + z)*((x + z + 1)*(-y + z + 1) - 4*z)
131  + (x + z)*(y - z)*(x - y + 2*z - 2))
132  - (x - y);
133 
134  out[2][0][0] = 0.5*(y - z)*(y - z + 1)*(2*x + 2*z - 1);
135  out[2][0][1] = 0.5*(x + z)*(2*y - 2*z + 1)*(x + z - 1);
136  out[2][0][2] = 0.5*(out[2][0][0] + out[2][0][1])
137  + 0.25*((y - x - 2*z)*(y - z + 1)*(x + z - 1)
138  + (x + z)*(y - z)*(y - x - 2*z + 2));
139 
140  out[3][0][0] = 0.5*(y - z)*(2*x + 2*z + 1)*(y - z + 1);
141  out[3][0][1] = 0.5*(2*y - 2*z + 1)*(x + z)*(x + z + 1);
142  out[3][0][2] = 0.5*(out[3][0][0] + out[3][0][1])
143  + 0.25*((y - x - 2*z)*(y - z + 1)*(x + z + 1)
144  + (y - z)*(x + z)*(y - x - 2*z));
145 
146  out[4][0][0] = 0;
147  out[4][0][1] = 0;
148  out[4][0][2] = 4*z - 1;
149 
150  // lower edges
151  out[5][0][0] = -((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
152  + (y - z + 1)*(x + z - 1)*((y - 1) + z));
153  out[5][0][1] = -((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
154  + (y - z + 1)*(x + z - 1)*((x + 1) - z));
155  out[5][0][2] = 0.5*(out[5][0][0] + out[5][0][1])
156  - 0.5*((-x + y - 2*z + 2)*((y - 1)*(x + 1) + z*(x - y + z + 1))
157  + (y - z + 1)*(x + z - 1)*(x - y + 2*z + 1));
158 
159  out[6][0][0] = -(y - z + 1)*(2*x + z + 1)*(y - 1);
160  out[6][0][1] = -(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1)
161  + (y - z + 1)*((x + z + 1)*x + 2*z));
162  out[6][0][2] = 0.5*(out[6][0][0] + out[6][0][1])
163  - 0.5*(-(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1))
164  + (y - z + 1)*(((y - 1)*x - 1) + 2*y + 1));
165 
166  out[7][0][0] = -(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1)
167  + (x + z - 1)*((y - z - 1)*y + 2*z));
168  out[7][0][1] = -(x + z - 1)*(2*y - z - 1)*(x + 1);
169  out[7][0][2] = 0.5*(out[7][0][0] + out[7][0][1])
170  - 0.5*(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1)
171  + (x + z - 1)*((-(x + 1)*y - 1) + 2*x + 1));
172 
173  out[8][0][0] = -(y - z + 1)*(2*x + z)*y;
174  out[8][0][1] = -(2*y - z + 1)*(x + z - 1)*(x + 1);
175  out[8][0][2] = 0.5*(out[8][0][0] + out[8][0][1])
176  - 0.5*(-x + y - 2*z + 2)*(x + 1)*y;
177 
178  // upper edges
179  out[9][0][0] = 2*z*(y - z - 1);
180  out[9][0][1] = 2*z*(x + z - 1);
181  out[9][0][2] = 0.5*(out[9][0][0] + out[9][0][1])
182  + (x + z - 1)*(y - z - 1) + z*(-x + y - 2*z);
183 
184  out[10][0][0] = -2*z*(y - z - 1);
185  out[10][0][1] = -2*z*(x + z + 1);
186  out[10][0][2] = 0.5*(out[10][0][0] + out[10][0][1])
187  - ((x + z + 1)*(y - z - 1) + 4*z)
188  - z*(-x + y - 2*z + 2);
189 
190  out[11][0][0] = -2*z*(y - z + 1);
191  out[11][0][1] = -2*z*(x + z - 1);
192  out[11][0][2] = 0.5*(out[11][0][0] + out[11][0][1])
193  - (y - z + 1)*(x + z - 1) - z*(-x + y - 2*z + 2);
194 
195  out[12][0][0] = 2*z*(y - z + 1);
196  out[12][0][1] = 2*z*(x + z + 1);
197  out[12][0][2] = 0.5*(out[12][0][0] + out[12][0][1])
198  + (y - z + 1)*(x + z + 1) + z*(-x + y - 2*z);
199 
200  // base face
201  out[13][0][0] = 2*((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
202  + (y - z + 1)*(x + z - 1)*(y - 1 + z));
203  out[13][0][1] = 2*((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
204  + (y - z + 1)*(x + z - 1)*(x + 1 - z));
205  out[13][0][2] = 0.5*(out[13][0][0] + out[13][0][1])
206  + ((-x + y - 2*z + 2)*((y - 1)*(x + 1) + z*(x - y + z + 1))
207  + (y - z + 1)*(x + z - 1)*(x - y + 2*z + 1));
208  }
209  else
210  {
211  // vertices
212  out[0][0][0] = 0.5*(y + z)*(y + z - 1)*(2*x - 2*z - 1);
213  out[0][0][1] = 0.5*(2*y + 2*z - 1)*(x - z - 1)*(x - z);
214  out[0][0][2] = 0.5*(out[0][0][0] + out[0][0][1])
215  + 0.25*((2*y + 2*z - 1)*(x - z - 1)*(x - z)
216  + (y + z)*(y + z - 1)*(-2*x + 2*z + 1));
217 
218  out[1][0][0] = -0.5*(y + z)*(2*x - 2*z + 1)*(-y - z + 1);
219  out[1][0][1] = -0.5*(x - z)*(x - z + 1)*(-2*y - 2*z + 1);
220  out[1][0][2] = 0.5*(out[1][0][0] + out[1][0][1])
221  - 0.25*((x - y - 2*z)*(x - z + 1)*(-y - z + 1)
222  + (x - z)*(y + z)*(-x + y + 2*z - 2));
223 
224  out[2][0][0] = 0.5*((y + z)*((x - z - 1)*(y + z + 1) + 4*z)
225  + (x - z)*(y + z)*(y + z + 1) + 4*z);
226  out[2][0][1] = 0.5*((x - z)*((x - z - 1)*(y + z + 1) + 4*z)
227  + (x - z)*(y + z)*(x - z - 1) - 4*z);
228  out[2][0][2] = 0.5*(out[2][0][0] + out[2][0][1])
229  + 0.25*((x - y - 2*z)*((x - z - 1)*(y + z + 1) + 4*z)
230  + (x - z)*(y + z)*(x - y - 2*z + 2) + 4*(x - y));
231 
232  out[3][0][0] = 0.5*(y + z)*(2*x - 2*z + 1)*(y + z + 1);
233  out[3][0][1] = 0.5*(x - z)*(x - z + 1)*(2*y + 2*z + 1);
234  out[3][0][2] = 0.5*(out[3][0][0] + out[3][0][1])
235  + 0.25*((x - y - 2*z)*(x - z + 1)*(y + z + 1)
236  + (y + z)*(x - z)*(x - y - 2*z));
237 
238  out[4][0][0] = 0;
239  out[4][0][1] = 0;
240  out[4][0][2] = 4*z - 1;
241 
242  // lower edges
243  out[5][0][0] = -(y + z - 1)*(2*x - z - 1)*(y + 1);
244  out[5][0][1] = -(((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1)
245  + (y + z - 1)*((x - z - 1)*x + 2*z));
246  out[5][0][2] = 0.5*(out[5][0][0] + out[5][0][1])
247  - 0.5*((((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1))
248  + (y + z - 1)*((-(y + 1)*x - 1) + 2*y + 1));
249 
250  out[6][0][0] = -(2*x - z + 1)*(y + z - 1)*(y + 1);
251  out[6][0][1] = -(x - z + 1)*(2*y + z)*x;
252  out[6][0][2] = 0.5*(out[6][0][0] + out[6][0][1])
253  - 0.5*(x - y - 2*z + 2)*(y + 1)*x;
254 
255  out[7][0][0] = -(2*x - z)*(y + z - 1)*y;
256  out[7][0][1] = -(x - z + 1)*(2*y + z - 1)*(x - 1);
257  out[7][0][2] = 0.5*(out[7][0][0] + out[7][0][1])
258  - 0.5*(x - y - 2*z + 2)*(x - 1)*y;
259 
260  out[8][0][0] = -(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1)
261  + (x - z + 1)*((y + z + 1)*y + 2*z));
262  out[8][0][1] = -(x - z + 1)*(2*y + z + 1)*(x - 1);
263  out[8][0][2] = 0.5*(out[8][0][0] + out[8][0][1])
264  - 0.5*(-(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1))
265  + (x - z + 1)*(((x - 1)*y - 1) + 2*x + 1));
266 
267  // upper edges
268  out[9][0][0] = 2*z*(y + z - 1);
269  out[9][0][1] = 2*z*(x - z - 1);
270  out[9][0][2] = 0.5*(out[9][0][0] + out[9][0][1])
271  + (y + z - 1)*(x - z - 1) + z*(x - y - 2*z);
272 
273  out[10][0][0] = -2*z*(y + z - 1);
274  out[10][0][1] = -2*z*(x - z + 1);
275  out[10][0][2] = 0.5*(out[10][0][0] + out[10][0][1])
276  - (x - z + 1)*(y + z - 1) - z*(x - y - 2*z + 2);
277 
278  out[11][0][0] = -2*z*(y + z + 1);
279  out[11][0][1] = -2*z*(x - z - 1);
280  out[11][0][2] = 0.5*(out[11][0][0] + out[11][0][1])
281  - ((y + z + 1)*(x - z - 1) + 4*z) - z*(x - y - 2*z + 2);
282 
283  out[12][0][0] = 2*z*(y + z + 1);
284  out[12][0][1] = 2*z*(x - z + 1);
285  out[12][0][2] = 0.5*(out[12][0][0] + out[12][0][1])
286  + (x - z + 1)*(y + z + 1) + z*(x - y - 2*z);
287 
288  // base face
289  out[13][0][0] = 2*((y + z - 1)*((y + 1)*(x - 1) - z*(x - y - z - 1))
290  + (x - z + 1)*(y + z - 1)*(y + 1 - z));
291  out[13][0][1] = 2*((x - z + 1)*((y + 1)*(x - 1) - z*(x - y - z - 1))
292  + (x - z + 1)*(y + z - 1)*(x - 1 + z));
293  out[13][0][2] = 0.5*(out[13][0][0] + out[13][0][1])
294  + (x - y - 2*z + 2)*((y + 1)*(x - 1) - z*(x - y - z - 1))
295  + (x - z + 1)*(y + z - 1)*(-(x - y - 2*z - 1));
296  }
297  }
298 
300  unsigned int order () const
301  {
302  return 2;
303  }
304  };
305 }
306 #endif