dune-localfunctions
2.2.0
|
00001 #ifndef DUNE_PK2DLOCALBASIS_HH 00002 #define DUNE_PK2DLOCALBASIS_HH 00003 00004 #include <dune/common/fmatrix.hh> 00005 00006 #include <dune/localfunctions/common/localbasis.hh> 00007 00008 namespace Dune 00009 { 00022 template<class D, class R, unsigned int k> 00023 class Pk2DLocalBasis 00024 { 00025 public: 00026 00028 enum {N = (k+1)*(k+2)/2}; 00029 00033 enum {O = k}; 00034 00035 typedef LocalBasisTraits<D,2,Dune::FieldVector<D,2>,R,1,Dune::FieldVector<R,1>, 00036 Dune::FieldMatrix<R,1,2> > Traits; 00037 00039 Pk2DLocalBasis () 00040 { 00041 for (unsigned int i=0; i<=k; i++) 00042 pos[i] = (1.0*i)/std::max(k,(unsigned int)1); 00043 } 00044 00046 unsigned int size () const 00047 { 00048 return N; 00049 } 00050 00052 inline void evaluateFunction (const typename Traits::DomainType& x, 00053 std::vector<typename Traits::RangeType>& out) const 00054 { 00055 out.resize(N); 00056 // specialization for k==0, not clear whether that is needed 00057 if (k==0) { 00058 out[0] = 1; 00059 return; 00060 } 00061 00062 int n=0; 00063 for (unsigned int j=0; j<=k; j++) 00064 for (unsigned int i=0; i<=k-j; i++) 00065 { 00066 out[n] = 1.0; 00067 for (unsigned int alpha=0; alpha<i; alpha++) 00068 out[n] *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]); 00069 for (unsigned int beta=0; beta<j; beta++) 00070 out[n] *= (x[1]-pos[beta])/(pos[j]-pos[beta]); 00071 for (unsigned int gamma=i+j+1; gamma<=k; gamma++) 00072 out[n] *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]); 00073 n++; 00074 } 00075 } 00076 00078 inline void 00079 evaluateJacobian (const typename Traits::DomainType& x, // position 00080 std::vector<typename Traits::JacobianType>& out) const // return value 00081 { 00082 out.resize(N); 00083 00084 // specialization for k==0, not clear whether that is needed 00085 if (k==0) { 00086 out[0][0][0] = 0; out[0][0][1] = 0; 00087 return; 00088 } 00089 00090 int n=0; 00091 for (unsigned int j=0; j<=k; j++) 00092 for (unsigned int i=0; i<=k-j; i++) 00093 { 00094 // x_0 derivative 00095 out[n][0][0] = 0.0; 00096 R factor=1.0; 00097 for (unsigned int beta=0; beta<j; beta++) 00098 factor *= (x[1]-pos[beta])/(pos[j]-pos[beta]); 00099 for (unsigned int a=0; a<i; a++) 00100 { 00101 R product=factor; 00102 for (unsigned int alpha=0; alpha<i; alpha++) 00103 if (alpha==a) 00104 product *= 1.0/(pos[i]-pos[alpha]); 00105 else 00106 product *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]); 00107 for (unsigned int gamma=i+j+1; gamma<=k; gamma++) 00108 product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]); 00109 out[n][0][0] += product; 00110 } 00111 for (unsigned int c=i+j+1; c<=k; c++) 00112 { 00113 R product=factor; 00114 for (unsigned int alpha=0; alpha<i; alpha++) 00115 product *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]); 00116 for (unsigned int gamma=i+j+1; gamma<=k; gamma++) 00117 if (gamma==c) 00118 product *= -1.0/(pos[gamma]-pos[i]-pos[j]); 00119 else 00120 product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]); 00121 out[n][0][0] += product; 00122 } 00123 00124 // x_1 derivative 00125 out[n][0][1] = 0.0; 00126 factor = 1.0; 00127 for (unsigned int alpha=0; alpha<i; alpha++) 00128 factor *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]); 00129 for (unsigned int b=0; b<j; b++) 00130 { 00131 R product=factor; 00132 for (unsigned int beta=0; beta<j; beta++) 00133 if (beta==b) 00134 product *= 1.0/(pos[j]-pos[beta]); 00135 else 00136 product *= (x[1]-pos[beta])/(pos[j]-pos[beta]); 00137 for (unsigned int gamma=i+j+1; gamma<=k; gamma++) 00138 product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]); 00139 out[n][0][1] += product; 00140 } 00141 for (unsigned int c=i+j+1; c<=k; c++) 00142 { 00143 R product=factor; 00144 for (unsigned int beta=0; beta<j; beta++) 00145 product *= (x[1]-pos[beta])/(pos[j]-pos[beta]); 00146 for (unsigned int gamma=i+j+1; gamma<=k; gamma++) 00147 if (gamma==c) 00148 product *= -1.0/(pos[gamma]-pos[i]-pos[j]); 00149 else 00150 product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]); 00151 out[n][0][1] += product; 00152 } 00153 00154 n++; 00155 } 00156 00157 } 00158 00160 unsigned int order () const 00161 { 00162 return k; 00163 } 00164 00165 private: 00166 R pos[k+1]; // positions on the interval 00167 }; 00168 00169 } 00170 #endif