1 #ifndef DUNE_Pk1DLOCALBASIS_HH
2 #define DUNE_Pk1DLOCALBASIS_HH
4 #include <dune/common/fmatrix.hh>
22 template<
class D,
class R,
unsigned int k>
35 Dune::FieldVector<D,1>,
38 Dune::FieldVector<R,1>,
39 Dune::FieldMatrix<R,1,1>
45 for (
unsigned int i=0; i<=k; i++)
46 pos[i] = (1.0*i)/std::max(k,(
unsigned int)1);
57 std::vector<typename Traits::RangeType>& out)
const
61 for (
unsigned int i=0; i<
N; i++)
64 for (
unsigned int alpha=0; alpha<i; alpha++)
65 out[i] *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
66 for (
unsigned int gamma=i+1; gamma<=k; gamma++)
67 out[i] *= (x[0]-pos[gamma])/(pos[i]-pos[gamma]);
74 std::vector<typename Traits::JacobianType>& out)
const
78 for (
unsigned int i=0; i<=k; i++) {
83 for (
unsigned int a=0; a<i; a++)
86 for (
unsigned int alpha=0; alpha<i; alpha++)
87 product *= (alpha==a) ? 1.0/(pos[i]-pos[alpha])
88 : (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
89 for (
unsigned int gamma=i+1; gamma<=k; gamma++)
90 product *= (pos[gamma]-x[0])/(pos[gamma]-pos[i]);
91 out[i][0][0] += product;
93 for (
unsigned int c=i+1; c<=k; c++)
96 for (
unsigned int alpha=0; alpha<i; alpha++)
97 product *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
98 for (
unsigned int gamma=i+1; gamma<=k; gamma++)
99 product *= (gamma==c) ? -1.0/(pos[gamma]-pos[i])
100 : (pos[gamma]-x[0])/(pos[gamma]-pos[i]);
101 out[i][0][0] += product;