1 #ifndef DUNE_PK2DLOCALBASIS_HH
2 #define DUNE_PK2DLOCALBASIS_HH
4 #include <dune/common/fmatrix.hh>
22 template<
class D,
class R,
unsigned int k>
28 enum {
N = (k+1)*(k+2)/2};
41 for (
unsigned int i=0; i<=k; i++)
42 pos[i] = (1.0*i)/std::max(k,(
unsigned int)1);
53 std::vector<typename Traits::RangeType>& out)
const
63 for (
unsigned int j=0; j<=k; j++)
64 for (
unsigned int i=0; i<=k-j; i++)
67 for (
unsigned int alpha=0; alpha<i; alpha++)
68 out[n] *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
69 for (
unsigned int beta=0; beta<j; beta++)
70 out[n] *= (x[1]-pos[beta])/(pos[j]-pos[beta]);
71 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
72 out[n] *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
80 std::vector<typename Traits::JacobianType>& out)
const
86 out[0][0][0] = 0; out[0][0][1] = 0;
91 for (
unsigned int j=0; j<=k; j++)
92 for (
unsigned int i=0; i<=k-j; i++)
97 for (
unsigned int beta=0; beta<j; beta++)
98 factor *= (x[1]-pos[beta])/(pos[j]-pos[beta]);
99 for (
unsigned int a=0; a<i; a++)
102 for (
unsigned int alpha=0; alpha<i; alpha++)
104 product *= 1.0/(pos[i]-pos[alpha]);
106 product *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
107 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
108 product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
109 out[n][0][0] += product;
111 for (
unsigned int c=i+j+1; c<=k; c++)
114 for (
unsigned int alpha=0; alpha<i; alpha++)
115 product *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
116 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
118 product *= -1.0/(pos[gamma]-pos[i]-pos[j]);
120 product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
121 out[n][0][0] += product;
127 for (
unsigned int alpha=0; alpha<i; alpha++)
128 factor *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
129 for (
unsigned int b=0; b<j; b++)
132 for (
unsigned int beta=0; beta<j; beta++)
134 product *= 1.0/(pos[j]-pos[beta]);
136 product *= (x[1]-pos[beta])/(pos[j]-pos[beta]);
137 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
138 product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
139 out[n][0][1] += product;
141 for (
unsigned int c=i+j+1; c<=k; c++)
144 for (
unsigned int beta=0; beta<j; beta++)
145 product *= (x[1]-pos[beta])/(pos[j]-pos[beta]);
146 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
148 product *= -1.0/(pos[gamma]-pos[i]-pos[j]);
150 product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
151 out[n][0][1] += product;