4 #ifndef DUNE_LOCALFUNCTIONS_QKLOCALBASIS_HH
5 #define DUNE_LOCALFUNCTIONS_QKLOCALBASIS_HH
7 #include <dune/common/fvector.hh>
8 #include <dune/common/fmatrix.hh>
9 #include <dune/common/power.hh>
11 #include <dune/geometry/type.hh>
31 template<
class D,
class R,
int k,
int d>
34 enum { n = StaticPower<k+1,d>::power };
37 static R p (
int i, D x)
40 for (
int j=0; j<=k; j++)
41 if (j!=i) result *= (k*x-j)/(i-j);
46 static R dp (
int i, D x)
50 for (
int j=0; j<=k; j++)
53 R prod( (k*1.0)/(i-j) );
54 for (
int l=0; l<=k; l++)
56 prod *= (k*x-l)/(i-l);
63 static Dune::FieldVector<int,d> multiindex (
int i)
65 Dune::FieldVector<int,d> alpha;
66 for (
int j=0; j<d; j++)
75 typedef LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d> >
Traits;
80 return StaticPower<k+1,d>::power;
85 std::vector<typename Traits::RangeType>& out)
const
88 for (
size_t i=0; i<
size(); i++)
91 Dune::FieldVector<int,d> alpha(multiindex(i));
97 for (
int j=0; j<d; j++)
98 out[i] *= p(alpha[j],in[j]);
108 std::vector<typename Traits::JacobianType>& out)
const
113 for (
size_t i=0; i<
size(); i++)
116 Dune::FieldVector<int,d> alpha(multiindex(i));
119 for (
int j=0; j<d; j++)
123 out[i][0][j] = dp(alpha[j],in[j]);
126 for (
int l=0; l<d; l++)
128 out[i][0][j] *= p(alpha[l],in[l]);