pk1dlocalbasis.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_Pk1DLOCALBASIS_HH
#define DUNE_Pk1DLOCALBASIS_HH
 
#include <dune/common/fmatrix.hh>
 
#include <dune/localfunctions/common/localbasis.hh>
 
namespace Dune
{
  template<class D, class R, unsigned int k>
  class Pk1DLocalBasis
  {
  public:
 
    enum {N = k+1};
 
    enum {O = k};
 
    typedef LocalBasisTraits<D,
        1,
        Dune::FieldVector<D,1>,
        R,
        1,
        Dune::FieldVector<R,1>,
        Dune::FieldMatrix<R,1,1>
        > Traits;
 
    Pk1DLocalBasis ()
    {
      for (unsigned int i=0; i<=k; i++)
        pos[i] = (1.0*i)/std::max(k,(unsigned int)1);
    }
 
    unsigned int size () const
    {
      return N;
    }
 
    inline void evaluateFunction (const typename Traits::DomainType& x,
                                  std::vector<typename Traits::RangeType>& out) const
    {
      out.resize(N);
 
      for (unsigned int i=0; i<N; i++)
      {
        out[i] = 1.0;
        for (unsigned int alpha=0; alpha<i; alpha++)
          out[i] *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
        for (unsigned int gamma=i+1; gamma<=k; gamma++)
          out[i] *= (x[0]-pos[gamma])/(pos[i]-pos[gamma]);
      }
    }
 
    inline void
    evaluateJacobian (const typename Traits::DomainType& x,             // position
                      std::vector<typename Traits::JacobianType>& out) const          // return value
    {
      out.resize(N);
 
      for (unsigned int i=0; i<=k; i++) {
 
        // x_0 derivative
        out[i][0][0] = 0.0;
        R factor=1.0;
        for (unsigned int a=0; a<i; a++)
        {
          R product=factor;
          for (unsigned int alpha=0; alpha<i; alpha++)
            product *= (alpha==a) ? 1.0/(pos[i]-pos[alpha])
                       : (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
          for (unsigned int gamma=i+1; gamma<=k; gamma++)
            product *= (pos[gamma]-x[0])/(pos[gamma]-pos[i]);
          out[i][0][0] += product;
        }
        for (unsigned int c=i+1; c<=k; c++)
        {
          R product=factor;
          for (unsigned int alpha=0; alpha<i; alpha++)
            product *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
          for (unsigned int gamma=i+1; gamma<=k; gamma++)
            product *= (gamma==c) ? -1.0/(pos[gamma]-pos[i])
                       : (pos[gamma]-x[0])/(pos[gamma]-pos[i]);
          out[i][0][0] += product;
        }
      }
 
    }
 
    void partial(const std::array<unsigned int,1>& order,
                 const typename Traits::DomainType& in,
                 std::vector<typename Traits::RangeType>& out) const
    {
      switch (order[0])
      {
        case 0:
          evaluateFunction(in,out);
          break;
        default:
          DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
      }
    }
    unsigned int order () const
    {
      return k;
    }
 
  private:
    R pos[k+1];     // positions on the interval
  };
 
}
#endif