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