dune-localfunctions  2.2.0
q2localbasis.hh
Go to the documentation of this file.
00001 // -*- tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*-
00002 // vi: set et ts=2 sw=2 sts=2:
00003 #ifndef DUNE_Q2_LOCALBASIS_HH
00004 #define DUNE_Q2_LOCALBASIS_HH
00005 
00006 #include <dune/common/fmatrix.hh>
00007 
00008 #include <dune/localfunctions/common/localbasis.hh>
00009 
00010 namespace Dune 
00011 {
00023 template<class D, class R, int dim>
00024 class Q2LocalBasis
00025 {
00026 public:
00027   typedef LocalBasisTraits<D,dim,Dune::FieldVector<D,dim>,R,1,Dune::FieldVector<R,1>,
00028       Dune::FieldMatrix<R,1,dim> > Traits;
00029 
00031   unsigned int size () const
00032   {
00033       int size = 1;
00034       for (int i=0; i<dim; i++)
00035           size *= 3;
00036       return size;
00037   }
00038 
00040   inline void evaluateFunction (const typename Traits::DomainType& in,
00041       std::vector<typename Traits::RangeType>& out) const
00042   {
00043     out.resize(size());
00044 
00045     // Evaluate the Lagrange functions
00046     array<array<R,3>, dim> X;
00047     
00048     for (size_t i=0; i<dim; i++) {
00049             X[i][0] =  R(2)*in[i]*in[i] - R(3)*in[i]+R(1);
00050             X[i][1] = -R(4)*in[i]*in[i] + R(4)*in[i];
00051             X[i][2] =  R(2)*in[i]*in[i] -   in[i];
00052     }
00053         
00054 
00055     // Compute function values: they are products of 1d Lagrange function values
00056     for (size_t i=0; i<out.size(); i++) {
00057             
00058         out[i] = 1;
00059                 
00060         // Construct the i-th Lagrange point
00061         size_t ternary = i;
00062         for (int j=0; j<dim; j++) {
00063                     
00064             int digit = ternary%3;
00065             ternary /= 3;
00066                     
00067             // Multiply the 1d Lagrange shape functions together
00068             out[i] *= X[j][digit];
00069                     
00070         }
00071                     
00072     }
00073 
00074   }
00075 
00077   inline void
00078   evaluateJacobian (const typename Traits::DomainType& in,         // position
00079       std::vector<typename Traits::JacobianType>& out) const      // return value
00080   {
00081     out.resize(size());
00082     
00083     // Evaluate the 1d Lagrange functions and their derivatives
00084     array<array<R,3>, dim> X, DX;
00085     
00086     for (size_t i=0; i<dim; i++) {
00087             X[i][0] =  R(2)*in[i]*in[i] - R(3)*in[i]+R(1);
00088             X[i][1] = -R(4)*in[i]*in[i] + R(4)*in[i];
00089             X[i][2] =  R(2)*in[i]*in[i] -   in[i];
00090             
00091             DX[i][0] =  R(4)*in[i] - R(3);
00092             DX[i][1] = -R(8)*in[i] + R(4);
00093             DX[i][2] =  R(4)*in[i] - R(1);
00094     }
00095         
00096     
00097     // Compute the derivatives by deriving the products of 1d Lagrange functions
00098     for (size_t i=0; i<out.size(); i++) {
00099             
00100         // Computing the j-th partial derivative
00101         for (int j=0; j<dim; j++) {
00102             
00103             out[i][0][j] = 1;
00104                 
00105             // Loop over the 'dim' terms in the product rule
00106             size_t ternary = i;
00107             for (int k=0; k<dim; k++) {
00108                     
00109                 int digit = ternary%3;
00110                 ternary /= 3;
00111                 
00112                 out[i][0][j] *= (k==j) ? DX[k][digit] : X[k][digit];
00113                     
00114             }
00115             
00116         }
00117                     
00118     }
00119 
00120 
00121 
00122   }
00123 
00125   unsigned int order () const
00126   {
00127     return 2;
00128   }
00129 };
00130 }
00131 #endif