dune-localfunctions  2.2.0
pk1dlocalbasis.hh
Go to the documentation of this file.
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