1 #ifndef DUNE_PK3DLOCALBASIS_HH
2 #define DUNE_PK3DLOCALBASIS_HH
4 #include <dune/common/fmatrix.hh>
22 template<
class D,
class R,
unsigned int k>
26 enum {
N = (k+1)*(k+2)*(k+3)/6};
43 std::vector<typename Traits::RangeType>& out)
const
51 for (i[2] = 0; i[2] <= k; ++i[2])
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])
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])
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];
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];
80 std::vector<typename Traits::JacobianType>& out)
const
88 for (i[2] = 0; i[2] <= k; ++i[2])
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])
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])
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];
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)
113 for (
unsigned int l = 0; l < i[3]; ++l)
120 for (
unsigned int j = 0; j < i[3]; ++j)
121 factor[3] *= kx3 - j;
122 for (
unsigned int m = 0; m < 3; ++m)
125 for (
unsigned int j = 0; j < i[m]; ++j)
128 for (
unsigned int p = 0; p < 3; ++p)
131 for (
unsigned int l = 0; l < i[p]; ++l)
133 prod *= R(k) / (i[p]-l);
135 prod *= (kx[p]-l) / (i[p]-l);
139 out[n][0][m] += prod;
157 template<
class D,
class R>
176 std::vector<typename Traits::RangeType>& out)
const
185 std::vector<typename Traits::JacobianType>& out)
const
194 template<
typename E,
typename F,
typename C>
195 void interpolate (
const E& e,
const F& f, std::vector<C>& out)
const