Dune Core Modules (2.6.0)

p23dlocalbasis.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_P2_3DLOCALBASIS_HH
4 #define DUNE_P2_3DLOCALBASIS_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 {
24  template<class D, class R>
25  class P23DLocalBasis
26  {
27  public:
29  typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,R,1,Dune::FieldVector<R,1>,
30  Dune::FieldMatrix<R,1,3> > Traits;
31 
33  unsigned int size () const
34  {
35  return 10;
36  }
37 
39  inline void evaluateFunction (const typename Traits::DomainType& in,
40  std::vector<typename Traits::RangeType>& out) const
41  {
42  out.resize(10);
43 
44  int coeff;
45  R a[2], b[3], c[3];
46 
47  // case 0:
48  coeff=2;
49  a[0]=1.0;
50  a[1]=0.5;
51  b[0]=-1.0;
52  b[1]=-1.0;
53  b[2]=-1.0;
54  c[0]=-1.0;
55  c[1]=-1.0;
56  c[2]=-1.0;
57 
58  out[0] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
59 
60  // case 1:
61  coeff=2;
62  a[0]=0.0;
63  a[1]=-0.5;
64  b[0]=1.0;
65  b[1]=0.0;
66  b[2]=0.0;
67  c[0]=1.0;
68  c[1]=0.0;
69  c[2]=0.0;
70 
71  out[1] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
72 
73  // case 2:
74  coeff=2;
75  a[0]=0.0;
76  a[1]=-0.5;
77  b[0]=0.0;
78  b[1]=1.0;
79  b[2]=0.0;
80  c[0]=0.0;
81  c[1]=1.0;
82  c[2]=0.0;
83 
84  out[2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
85 
86  // case 3:
87  coeff=2;
88  a[0]=0.0;
89  a[1]=-0.5;
90  b[0]=0.0;
91  b[1]=0.0;
92  b[2]=1.0;
93  c[0]=0.0;
94  c[1]=0.0;
95  c[2]=1.0;
96 
97  out[3] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
98 
99  // case 4:
100  coeff=4;
101  a[0]=0.0;
102  a[1]=1.0;
103  b[0]=1.0;
104  b[1]=0.0;
105  b[2]=0.0;
106  c[0]=-1.0;
107  c[1]=-1.0;
108  c[2]=-1.0;
109 
110  out[4] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
111 
112  // case 5:
113  coeff=4;
114  a[0]=0.0;
115  a[1]=0.0;
116  b[0]=1.0;
117  b[1]=0.0;
118  b[2]=0.0;
119  c[0]=0.0;
120  c[1]=1.0;
121  c[2]=0.0;
122 
123  out[5] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
124 
125  // case 6:
126  coeff=4;
127  a[0]=0.0;
128  a[1]=1.0;
129  b[0]=0.0;
130  b[1]=1.0;
131  b[2]=0.0;
132  c[0]=-1.0;
133  c[1]=-1.0;
134  c[2]=-1.0;
135 
136  out[6] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
137 
138  // case 7:
139  coeff=4;
140  a[0]=0.0;
141  a[1]=1.0;
142  b[0]=0.0;
143  b[1]=0.0;
144  b[2]=1.0;
145  c[0]=-1.0;
146  c[1]=-1.0;
147  c[2]=-1.0;
148 
149  out[7] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
150 
151  // case 8:
152  coeff=4;
153  a[0]=0.0;
154  a[1]=0.0;
155  b[0]=1.0;
156  b[1]=0.0;
157  b[2]=0.0;
158  c[0]=0.0;
159  c[1]=0.0;
160  c[2]=1.0;
161 
162  out[8] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
163 
164  // case 9:
165  coeff=4;
166  a[0]=0.0;
167  a[1]=0.0;
168  b[0]=0.0;
169  b[1]=1.0;
170  b[2]=0.0;
171  c[0]=0.0;
172  c[1]=0.0;
173  c[2]=1.0;
174 
175  out[9] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
176 
177  }
178 
180  inline void
181  evaluateJacobian (const typename Traits::DomainType& in, // position
182  std::vector<typename Traits::JacobianType>& out) const // return value
183  {
184  out.resize(10);
185 
186  R aa[3][3], bb[3][3];
187  // case 0:
188  //x derivative
189  aa[0][0]=-3.0;
190  bb[0][0]=4.0;
191  bb[1][0]=4.0;
192  bb[2][0]=4.0;
193  //y derivative
194  aa[0][1]=-3.0;
195  bb[0][1]=4.0;
196  bb[1][1]=4.0;
197  bb[2][1]=4.0;
198  // z derivative
199  aa[0][2]=-3.0;
200  bb[0][2]=4.0;
201  bb[1][2]=4.0;
202  bb[2][2]=4.0;
203 
204  out[0][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
205  out[0][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
206  out[0][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
207 
208  //case 1:
209  //x derivative
210  aa[0][0]=-1.0;
211  bb[0][0]=4.0;
212  bb[1][0]=0.0;
213  bb[2][0]=0.0;
214  //y derivative
215  aa[0][1]=0.0;
216  bb[0][1]=0.0;
217  bb[1][1]=0.0;
218  bb[2][1]=0.0;
219  // z derivative
220  aa[0][2]=0.0;
221  bb[0][2]=0.0;
222  bb[1][2]=0.0;
223  bb[2][2]=0.0;
224 
225  out[1][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
226  out[1][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
227  out[1][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
228 
229  // case 2:
230  //x derivative
231  aa[0][0]=0.0;
232  bb[0][0]=0.0;
233  bb[1][0]=0.0;
234  bb[2][0]=0.0;
235  //y derivative
236  aa[0][1]=-1.0;
237  bb[0][1]=0.0;
238  bb[1][1]=4.0;
239  bb[2][1]=0.0;
240  // z derivative
241  aa[0][2]=0.0;
242  bb[0][2]=0.0;
243  bb[1][2]=0.0;
244  bb[2][2]=0.0;
245 
246  out[2][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
247  out[2][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
248  out[2][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
249 
250  // case 3:
251  //x derivative
252  aa[0][0]=0.0;
253  bb[0][0]=0.0;
254  bb[1][0]=0.0;
255  bb[2][0]=0.0;
256  //y derivative
257  aa[0][1]=0.0;
258  bb[0][1]=0.0;
259  bb[1][1]=0.0;
260  bb[2][1]=0.0;
261  // z derivative
262  aa[0][2]=-1.0;
263  bb[0][2]=0.0;
264  bb[1][2]=0.0;
265  bb[2][2]=4.0;
266 
267  out[3][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
268  out[3][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
269  out[3][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
270 
271  // case 4:
272  //x derivative
273  aa[0][0]=4.0;
274  bb[0][0]=-8.0;
275  bb[1][0]=-4.0;
276  bb[2][0]=-4.0;
277  //y derivative
278  aa[0][1]=0.0;
279  bb[0][1]=-4.0;
280  bb[1][1]=0.0;
281  bb[2][1]=0.0;
282  // z derivative
283  aa[0][2]=0.0;
284  bb[0][2]=-4.0;
285  bb[1][2]=0.0;
286  bb[2][2]=0.0;
287 
288  out[4][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
289  out[4][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
290  out[4][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
291 
292  // case 5:
293  //x derivative
294  aa[0][0]=0.0;
295  bb[0][0]=0.0;
296  bb[1][0]=4.0;
297  bb[2][0]=0.0;
298  //y derivative
299  aa[0][1]=0.0;
300  bb[0][1]=4.0;
301  bb[1][1]=0.0;
302  bb[2][1]=0.0;
303  // z derivative
304  aa[0][2]=0.0;
305  bb[0][2]=0.0;
306  bb[1][2]=0.0;
307  bb[2][2]=0.0;
308 
309  out[5][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
310  out[5][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
311  out[5][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
312 
313  // case 6:
314  //x derivative
315  aa[0][0]=0.0;
316  bb[0][0]=0.0;
317  bb[1][0]=-4.0;
318  bb[2][0]=0.0;
319  //y derivative
320  aa[0][1]=4.0;
321  bb[0][1]=-4.0;
322  bb[1][1]=-8.0;
323  bb[2][1]=-4.0;
324  // z derivative
325  aa[0][2]=0.0;
326  bb[0][2]=0.0;
327  bb[1][2]=-4.0;
328  bb[2][2]=0.0;
329 
330  out[6][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
331  out[6][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
332  out[6][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
333 
334  // case 7:
335  //x derivative
336  aa[0][0]=0.0;
337  bb[0][0]=0.0;
338  bb[1][0]=0.0;
339  bb[2][0]=-4.0;
340  //y derivative
341  aa[0][1]=0.0;
342  bb[0][1]=0.0;
343  bb[1][1]=0.0;
344  bb[2][1]=-4.0;
345  // z derivative
346  aa[0][2]=4.0;
347  bb[0][2]=-4.0;
348  bb[1][2]=-4.0;
349  bb[2][2]=-8.0;
350 
351  out[7][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
352  out[7][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
353  out[7][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
354 
355  //case 8:
356  //x derivative
357  aa[0][0]=0.0;
358  bb[0][0]=0.0;
359  bb[1][0]=0.0;
360  bb[2][0]=4.0;
361  //y derivative
362  aa[0][1]=0.0;
363  bb[0][1]=0.0;
364  bb[1][1]=0.0;
365  bb[2][1]=0.0;
366  // z derivative
367  aa[0][2]=0.0;
368  bb[0][2]=4.0;
369  bb[1][2]=0.0;
370  bb[2][2]=0.0;
371 
372  out[8][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
373  out[8][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
374  out[8][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
375 
376  // case 9:
377  //x derivative
378  aa[0][0]=0.0;
379  bb[0][0]=0.0;
380  bb[1][0]=0.0;
381  bb[2][0]=0.0;
382  //y derivative
383  aa[0][1]=0.0;
384  bb[0][1]=0.0;
385  bb[1][1]=0.0;
386  bb[2][1]=4.0;
387  // z derivative
388  aa[0][2]=0.0;
389  bb[0][2]=0.0;
390  bb[1][2]=4.0;
391  bb[2][2]=0.0;
392 
393  out[9][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
394  out[9][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
395  out[9][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
396 
397  }
398 
404  void partial(const std::array<unsigned int,3>& order,
405  const typename Traits::DomainType& in,
406  std::vector<typename Traits::RangeType>& out) const
407  {
408  auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
409  if (totalOrder == 0) {
410  evaluateFunction(in, out);
411  } else if (totalOrder == 1) {
412  out.resize(size());
413 
414  auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
415  switch (direction) {
416  case 0:
417  out[0] =-3.0 + 4.0*(in[0] + in[1] + in[2]);
418  out[1] =-1.0 + 4.0*in[0];
419  out[2] = 0.0;
420  out[3] = 0.0;
421  out[4] = 4.0 - 4.0*(2.0*in[0] + in[1] + in[2]);
422  out[5] = 4.0*in[1];
423  out[6] =-4.0*in[1];
424  out[7] =-4.0*in[2];
425  out[8] = 4.0*in[2];
426  out[9] = 0.0;
427  break;
428  case 1:
429  out[0] =-3.0 + 4.0*(in[0] + in[1] + in[2]);
430  out[1] = 0.0;
431  out[2] =-1.0 + 4.0*in[1];
432  out[3] = 0.0;
433  out[4] =-4.0*in[0];
434  out[5] = 4.0*in[0];
435  out[6] = 4.0 - 4.0*(in[0] + 2.0*in[1] + in[2]);
436  out[7] =-4.0*in[2];
437  out[8] = 0.0;
438  out[9] = 4.0*in[2];
439  break;
440  case 2:
441  out[0] =-3.0 + 4.0*(in[0] + in[1] + in[2]);
442  out[1] = 0.0;
443  out[2] = 0.0;
444  out[3] =-1.0 + 4.0*in[2];
445  out[4] =-4.0*in[0];
446  out[5] = 0.0;
447  out[6] =-4.0*in[1];
448  out[7] = 4.0 - 4.0*(in[0] + in[1] + 2.0*in[2]);
449  out[8] = 4.0*in[0];
450  out[9] = 4.0*in[1];
451  break;
452  default:
453  DUNE_THROW(RangeError, "Component out of range.");
454  }
455  } else {
456  DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
457  }
458  }
459 
461  unsigned int order () const
462  {
463  return 2;
464  }
465  };
466 }
467 #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::P23DLocalBasis
Quadratic Lagrange shape functions on the tetrahedron.
Definition: p23dlocalbasis.hh:26
Dune::P23DLocalBasis::order
unsigned int order() const
Polynomial order of the shape functions.
Definition: p23dlocalbasis.hh:461
Dune::P23DLocalBasis::size
unsigned int size() const
number of shape functions
Definition: p23dlocalbasis.hh:33
Dune::P23DLocalBasis::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: p23dlocalbasis.hh:30
Dune::P23DLocalBasis::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 any order of all shape functions.
Definition: p23dlocalbasis.hh:404
Dune::P23DLocalBasis::evaluateJacobian
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: p23dlocalbasis.hh:181
Dune::P23DLocalBasis::evaluateFunction
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: p23dlocalbasis.hh:39
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 7, 22:32, 2024)