3 #ifndef DUNE_PK3DLOCALBASIS_HH
4 #define DUNE_PK3DLOCALBASIS_HH
6 #include <dune/common/fmatrix.hh>
24 template<
class D,
class R,
unsigned int k>
28 enum {
N = (k+1)*(k+2)*(k+3)/6};
45 std::vector<typename Traits::RangeType>& out)
const
53 for (i[2] = 0; i[2] <= k; ++i[2])
56 for (
unsigned int j = 0; j < i[2]; ++j)
57 factor[2] *= (kx[2]-j) / (i[2]-j);
58 for (i[1] = 0; i[1] <= k - i[2]; ++i[1])
61 for (
unsigned int j = 0; j < i[1]; ++j)
62 factor[1] *= (kx[1]-j) / (i[1]-j);
63 for (i[0] = 0; i[0] <= k - i[1] - i[2]; ++i[0])
66 for (
unsigned int j = 0; j < i[0]; ++j)
67 factor[0] *= (kx[0]-j) / (i[0]-j);
68 i[3] = k - i[0] - i[1] - i[2];
69 D kx3 = k - kx[0] - kx[1] - kx[2];
71 for (
unsigned int j = 0; j < i[3]; ++j)
72 factor[3] *= (kx3-j) / (i[3]-j);
73 out[n++] = factor[0] * factor[1] * factor[2] * factor[3];
82 std::vector<typename Traits::JacobianType>& out)
const
90 for (i[2] = 0; i[2] <= k; ++i[2])
93 for (
unsigned int j = 0; j < i[2]; ++j)
94 factor[2] *= (kx[2]-j) / (i[2]-j);
95 for (i[1] = 0; i[1] <= k - i[2]; ++i[1])
98 for (
unsigned int j = 0; j < i[1]; ++j)
99 factor[1] *= (kx[1]-j) / (i[1]-j);
100 for (i[0] = 0; i[0] <= k - i[1] - i[2]; ++i[0])
103 for (
unsigned int j = 0; j < i[0]; ++j)
104 factor[0] *= (kx[0]-j) / (i[0]-j);
105 i[3] = k - i[0] - i[1] - i[2];
106 D kx3 = k - kx[0] - kx[1] - kx[2];
109 for (
unsigned int j = 0; j < i[3]; ++j)
110 factor[3] /= i[3] - j;
111 R prod_all = factor[0] * factor[1] * factor[2] * factor[3];
112 for (
unsigned int j = 0; j < i[3]; ++j)
115 for (
unsigned int l = 0; l < i[3]; ++l)
122 for (
unsigned int j = 0; j < i[3]; ++j)
123 factor[3] *= kx3 - j;
124 for (
unsigned int m = 0; m < 3; ++m)
127 for (
unsigned int j = 0; j < i[m]; ++j)
130 for (
unsigned int p = 0; p < 3; ++p)
133 for (
unsigned int l = 0; l < i[p]; ++l)
135 prod *= R(k) / (i[p]-l);
137 prod *= (kx[p]-l) / (i[p]-l);
141 out[n][0][m] += prod;
159 template<
class D,
class R>
178 std::vector<typename Traits::RangeType>& out)
const
187 std::vector<typename Traits::JacobianType>& out)
const
196 template<
typename E,
typename F,
typename C>
197 void interpolate (
const E& e,
const F& f, std::vector<C>& out)
const