dune-localfunctions
2.2.0
|
00001 #ifndef DUNE_Pk1DLOCALBASIS_HH 00002 #define DUNE_Pk1DLOCALBASIS_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 Pk1DLocalBasis 00024 { 00025 public: 00026 00028 enum {N = k+1}; 00029 00031 enum {O = k}; 00032 00033 typedef LocalBasisTraits<D, 00034 1, 00035 Dune::FieldVector<D,1>, 00036 R, 00037 1, 00038 Dune::FieldVector<R,1>, 00039 Dune::FieldMatrix<R,1,1> 00040 > Traits; 00041 00043 Pk1DLocalBasis () 00044 { 00045 for (unsigned int i=0; i<=k; i++) 00046 pos[i] = (1.0*i)/std::max(k,(unsigned int)1); 00047 } 00048 00050 unsigned int size () const 00051 { 00052 return N; 00053 } 00054 00056 inline void evaluateFunction (const typename Traits::DomainType& x, 00057 std::vector<typename Traits::RangeType>& out) const 00058 { 00059 out.resize(N); 00060 00061 for (unsigned int i=0; i<N; i++) 00062 { 00063 out[i] = 1.0; 00064 for (unsigned int alpha=0; alpha<i; alpha++) 00065 out[i] *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]); 00066 for (unsigned int gamma=i+1; gamma<=k; gamma++) 00067 out[i] *= (x[0]-pos[gamma])/(pos[i]-pos[gamma]); 00068 } 00069 } 00070 00072 inline void 00073 evaluateJacobian (const typename Traits::DomainType& x, // position 00074 std::vector<typename Traits::JacobianType>& out) const // return value 00075 { 00076 out.resize(N); 00077 00078 for (unsigned int i=0; i<=k; i++) { 00079 00080 // x_0 derivative 00081 out[i][0][0] = 0.0; 00082 R factor=1.0; 00083 for (unsigned int a=0; a<i; a++) 00084 { 00085 R product=factor; 00086 for (unsigned int alpha=0; alpha<i; alpha++) 00087 product *= (alpha==a) ? 1.0/(pos[i]-pos[alpha]) 00088 : (x[0]-pos[alpha])/(pos[i]-pos[alpha]); 00089 for (unsigned int gamma=i+1; gamma<=k; gamma++) 00090 product *= (pos[gamma]-x[0])/(pos[gamma]-pos[i]); 00091 out[i][0][0] += product; 00092 } 00093 for (unsigned int c=i+1; c<=k; c++) 00094 { 00095 R product=factor; 00096 for (unsigned int alpha=0; alpha<i; alpha++) 00097 product *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]); 00098 for (unsigned int gamma=i+1; gamma<=k; gamma++) 00099 product *= (gamma==c) ? -1.0/(pos[gamma]-pos[i]) 00100 : (pos[gamma]-x[0])/(pos[gamma]-pos[i]); 00101 out[i][0][0] += product; 00102 } 00103 } 00104 00105 } 00106 00108 unsigned int order () const 00109 { 00110 return k; 00111 } 00112 00113 private: 00114 R pos[k+1]; // positions on the interval 00115 }; 00116 00117 } 00118 #endif