dune-localfunctions  2.2.1
pk3dlocalbasis.hh
Go to the documentation of this file.
1 #ifndef DUNE_PK3DLOCALBASIS_HH
2 #define DUNE_PK3DLOCALBASIS_HH
3 
4 #include <dune/common/fmatrix.hh>
5 
7 
8 namespace Dune
9 {
22 template<class D, class R, unsigned int k>
24 {
25 public:
26  enum {N = (k+1)*(k+2)*(k+3)/6};
27  enum {O = k};
28 
29  typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,R,1,Dune::FieldVector<R,1>,
30  Dune::FieldMatrix<R,1,3> > Traits;
31 
34 
36  unsigned int size () const
37  {
38  return N;
39  }
40 
42  inline void evaluateFunction (const typename Traits::DomainType& x,
43  std::vector<typename Traits::RangeType>& out) const
44  {
45  out.resize(N);
46  typename Traits::DomainType kx = x;
47  kx *= k;
48  unsigned int n = 0;
49  unsigned int i[4];
50  R factor[4];
51  for (i[2] = 0; i[2] <= k; ++i[2])
52  {
53  factor[2] = 1.0;
54  for (unsigned int j = 0; j < i[2]; ++j)
55  factor[2] *= (kx[2]-j) / (i[2]-j);
56  for (i[1] = 0; i[1] <= k - i[2]; ++i[1])
57  {
58  factor[1] = 1.0;
59  for (unsigned int j = 0; j < i[1]; ++j)
60  factor[1] *= (kx[1]-j) / (i[1]-j);
61  for (i[0] = 0; i[0] <= k - i[1] - i[2]; ++i[0])
62  {
63  factor[0] = 1.0;
64  for (unsigned int j = 0; j < i[0]; ++j)
65  factor[0] *= (kx[0]-j) / (i[0]-j);
66  i[3] = k - i[0] - i[1] - i[2];
67  D kx3 = k - kx[0] - kx[1] - kx[2];
68  factor[3] = 1.0;
69  for (unsigned int j = 0; j < i[3]; ++j)
70  factor[3] *= (kx3-j) / (i[3]-j);
71  out[n++] = factor[0] * factor[1] * factor[2] * factor[3];
72  }
73  }
74  }
75  }
76 
78  inline void
79  evaluateJacobian (const typename Traits::DomainType& x, // position
80  std::vector<typename Traits::JacobianType>& out) const // return value
81  {
82  out.resize(N);
83  typename Traits::DomainType kx = x;
84  kx *= k;
85  unsigned int n = 0;
86  unsigned int i[4];
87  R factor[4];
88  for (i[2] = 0; i[2] <= k; ++i[2])
89  {
90  factor[2] = 1.0;
91  for (unsigned int j = 0; j < i[2]; ++j)
92  factor[2] *= (kx[2]-j) / (i[2]-j);
93  for (i[1] = 0; i[1] <= k - i[2]; ++i[1])
94  {
95  factor[1] = 1.0;
96  for (unsigned int j = 0; j < i[1]; ++j)
97  factor[1] *= (kx[1]-j) / (i[1]-j);
98  for (i[0] = 0; i[0] <= k - i[1] - i[2]; ++i[0])
99  {
100  factor[0] = 1.0;
101  for (unsigned int j = 0; j < i[0]; ++j)
102  factor[0] *= (kx[0]-j) / (i[0]-j);
103  i[3] = k - i[0] - i[1] - i[2];
104  D kx3 = k - kx[0] - kx[1] - kx[2];
105  R sum3 = 0.0;
106  factor[3] = 1.0;
107  for (unsigned int j = 0; j < i[3]; ++j)
108  factor[3] /= i[3] - j;
109  R prod_all = factor[0] * factor[1] * factor[2] * factor[3];
110  for (unsigned int j = 0; j < i[3]; ++j)
111  {
112  R prod = prod_all;
113  for (unsigned int l = 0; l < i[3]; ++l)
114  if (j == l)
115  prod *= -R(k);
116  else
117  prod *= kx3 - l;
118  sum3 += prod;
119  }
120  for (unsigned int j = 0; j < i[3]; ++j)
121  factor[3] *= kx3 - j;
122  for (unsigned int m = 0; m < 3; ++m)
123  {
124  out[n][0][m] = sum3;
125  for (unsigned int j = 0; j < i[m]; ++j)
126  {
127  R prod = factor[3];
128  for (unsigned int p = 0; p < 3; ++p)
129  {
130  if (m == p)
131  for (unsigned int l = 0; l < i[p]; ++l)
132  if (j == l)
133  prod *= R(k) / (i[p]-l);
134  else
135  prod *= (kx[p]-l) / (i[p]-l);
136  else
137  prod *= factor[p];
138  }
139  out[n][0][m] += prod;
140  }
141  }
142  n++;
143  }
144  }
145  }
146  }
147 
149  unsigned int order () const
150  {
151  return k;
152  }
153 };
154 
155 
156 //Specialization for k=0
157 template<class D, class R>
158 class Pk3DLocalBasis<D,R,0>
159 {
160 public:
161  typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,R,1,Dune::FieldVector<R,1>,
162  Dune::FieldMatrix<R,1,3> > Traits;
163 
165  enum {N = 1};
166 
168  enum {O = 0};
169 
170  unsigned int size () const
171  {
172  return 1;
173  }
174 
175  inline void evaluateFunction (const typename Traits::DomainType& in,
176  std::vector<typename Traits::RangeType>& out) const
177  {
178  out.resize(1);
179  out[0] = 1;
180  }
181 
182  // evaluate derivative of a single component
183  inline void
184  evaluateJacobian (const typename Traits::DomainType& in, // position
185  std::vector<typename Traits::JacobianType>& out) const // return value
186  {
187  out.resize(1);
188  out[0][0][0] = 0;
189  out[0][0][1] = 0;
190  out[0][0][2] = 0;
191  }
192 
193  // local interpolation of a function
194  template<typename E, typename F, typename C>
195  void interpolate (const E& e, const F& f, std::vector<C>& out) const
196  {
197  typename Traits::DomainType x;
198  typename Traits::RangeType y;
199  x[0] = 1.0/4.0;
200  x[1] = 1.0/4.0;
201  x[2] = 1.0/4.0;
202  f.eval_local(e,x,y);
203  out[0] = y;
204  }
205 
206  unsigned int order () const
207  {
208  return 0;
209  }
210 };
211 }
212 #endif