qklocalbasis.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_LOCALFUNCTIONS_QKLOCALBASIS_HH
#define DUNE_LOCALFUNCTIONS_QKLOCALBASIS_HH
 
#include <numeric>
 
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/power.hh>
 
#include <dune/geometry/type.hh>
 
#include <dune/localfunctions/common/localbasis.hh>
#include <dune/localfunctions/common/localfiniteelementtraits.hh>
 
 
namespace Dune
{
  template<class D, class R, int k, int d>
  class QkLocalBasis
  {
    // ith Lagrange polynomial of degree k in one dimension
    static R p (int i, D x)
    {
      R result(1.0);
      for (int j=0; j<=k; j++)
        if (j!=i) result *= (k*x-j)/(i-j);
      return result;
    }
 
    // derivative of ith Lagrange polynomial of degree k in one dimension
    static R dp (int i, D x)
    {
      R result(0.0);
 
      for (int j=0; j<=k; j++)
        if (j!=i)
        {
          R prod( (k*1.0)/(i-j) );
          for (int l=0; l<=k; l++)
            if (l!=i && l!=j)
              prod *= (k*x-l)/(i-l);
          result += prod;
        }
      return result;
    }
 
    // Second derivative of j-th Lagrange polynomial of degree k in one dimension
    // Formula and notation taken from https://en.wikipedia.org/wiki/Lagrange_polynomial#Derivatives
    static R ddp (int j, D x)
    {
      R result(0.0);
 
      for (int i=0; i<=k; i++)
      {
        if (i==j)
          continue;
 
        R sum(0);
 
        for (int m=0; m<=k; m++)
        {
          if (m==i || m==j)
            continue;
 
          R prod( (k*1.0)/(j-m) );
          for (int l=0; l<=k; l++)
            if (l!=i && l!=j && l!=m)
              prod *= (k*x-l)/(j-l);
          sum += prod;
        }
 
        result += sum * ( (k*1.0)/(j-i) );
      }
 
      return result;
    }
 
    // Return i as a d-digit number in the (k+1)-nary system
    static Dune::FieldVector<int,d> multiindex (int i)
    {
      Dune::FieldVector<int,d> alpha;
      for (int j=0; j<d; j++)
      {
        alpha[j] = i % (k+1);
        i = i/(k+1);
      }
      return alpha;
    }
 
  public:
    typedef LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d> > Traits;
 
    unsigned int size () const
    {
      return StaticPower<k+1,d>::power;
    }
 
    inline void evaluateFunction (const typename Traits::DomainType& in,
                                  std::vector<typename Traits::RangeType>& out) const
    {
      out.resize(size());
      for (size_t i=0; i<size(); i++)
      {
        // convert index i to multiindex
        Dune::FieldVector<int,d> alpha(multiindex(i));
 
        // initialize product
        out[i] = 1.0;
 
        // dimension by dimension
        for (int j=0; j<d; j++)
          out[i] *= p(alpha[j],in[j]);
      }
    }
 
    inline void
    evaluateJacobian (const typename Traits::DomainType& in,
                      std::vector<typename Traits::JacobianType>& out) const
    {
      out.resize(size());
 
      // Loop over all shape functions
      for (size_t i=0; i<size(); i++)
      {
        // convert index i to multiindex
        Dune::FieldVector<int,d> alpha(multiindex(i));
 
        // Loop over all coordinate directions
        for (int j=0; j<d; j++)
        {
          // Initialize: the overall expression is a product
          // if j-th bit of i is set to -1, else 1
          out[i][0][j] = dp(alpha[j],in[j]);
 
          // rest of the product
          for (int l=0; l<d; l++)
            if (l!=j)
              out[i][0][j] *= p(alpha[l],in[l]);
        }
      }
    }
 
    inline void partial(const std::array<unsigned int,d>& order,
                        const typename Traits::DomainType& in,
                        std::vector<typename Traits::RangeType>& out) const
    {
      out.resize(size());
 
      // Loop over all shape functions
      for (size_t i=0; i<size(); i++)
      {
        // convert index i to multiindex
        Dune::FieldVector<int,d> alpha(multiindex(i));
 
        // Initialize: the overall expression is a product
        out[i][0] = 1.0;
 
        // rest of the product
        for (std::size_t l=0; l<d; l++)
        {
          switch (order[l])
          {
            case 0:
              out[i][0] *= p(alpha[l],in[l]);
              break;
            case 1:
              out[i][0] *= dp(alpha[l],in[l]);
              break;
            case 2:
              out[i][0] *= ddp(alpha[l],in[l]);
              break;
            default:
              DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
          }
        }
      }
    }
 
    unsigned int order () const
    {
      return k;
    }
  };
}
 
#endif