serendipitylocalbasis.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_LOCALFUNCTIONS_SERENDIPITYLOCALBASIS_HH
#define DUNE_LOCALFUNCTIONS_SERENDIPITYLOCALBASIS_HH
 
#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 SerendipityLocalBasis
  {
    typedef typename Dune::FieldVector<int,d> DuneVector;
 
 
    // ith Lagrange polynomial of degree k in one dimension
    static R p (int i,int deg, D x)
    {
      R result(1.0);
      for (int j=0; j<=deg; j++)
        if (j!=i) result *= (deg*x-j)/(i-j);
      return result;
    }
 
    // derivative of ith Lagrange polynomial of degree k in one dimension
    static R dp (int i, int deg, D x)
    {
      R result(0.0);
 
      for (int j=0; j<=deg; j++)
        if (j!=i)
        {
          R prod( (deg*1.0)/(i-j) );
          for (int l=0; l<=deg; l++)
            if (l!=i && l!=j)
              prod *= (deg*x-l)/(i-l);
          result += prod;
        }
      return result;
    }
 
    static void multiindex (DuneVector& alpha)
    {
            int i = 0;
            for (int j=0; j<d; j++){
                    i+=alpha[j]*pow(k+1,j);
            }
            i++;
            for (int j=0; j<d; j++)
            {
                    alpha[j] = i % (k+1);
                    i = i/(k+1);
            }
            unsigned int sum = 0;
            for (unsigned int j=0; j<d; j++){
          if(alpha[j] > 0 && alpha[j] < k){ sum++; }
            }
            if(sum > 1){
          multiindex(alpha);
            }
    }
 
    bool isEdge(DuneVector& alpha) const{
      for(int i = 0; i < d; i++){
        if(alpha[i]>0 && alpha[i] < k)
          return true;
      }
      return false;
    }
 
    double edgeVal(int alpha, const double& in, int der) const{
        //add contribution from alpha
        if(alpha > 0 && alpha < k){
          if(der == 0)
            return p(alpha,k,in);
          else
            return dp(alpha,k,in);
        }
        else{
          if(der == 0)
            return p(alpha/k,k-1,in);
          else
            return dp(alpha/k,k-1,in);
        }
    }
 
    double nodeVal(int alpha, const double& in, int der) const{
      // (k-1)-th order component
      if(k == 1 && der == 0)
        return p(alpha,k,in);
      if(k == 1 && der > 0)
        return dp(alpha,k, in); 
      if(der == 0)
        return p(alpha/k,k-1,in);
      else
        return dp(alpha/k,k-1,in);
    }
 
  public:
    typedef LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d>, 1> Traits;
 
    unsigned int size () const
    {
      if(d == 1){
        return k+1;
      }
      if(d == 2){
        return 2*(k+1)+2*(k-1);
      }
      if(d == 3){
        return 4*(k+1)+8*(k-1);
      }
    }
 
    inline void evaluateFunction (const typename Traits::DomainType& in,
                                  std::vector<typename Traits::RangeType>& out) const
    {
      DuneVector alpha(0);
      out.resize(size());
      for (size_t i=0; i<size(); i++)
      {
        // convert index i to multiindex
        if(i>0){ multiindex(alpha); }
 
        out[i] = 1.0;
        for(int j = 0; j < d; j++){
          if(isEdge(alpha))
            out[i] *= edgeVal(alpha[j], in[j], 0);
          else{ //isNode
            out[i] *= nodeVal(alpha[j], in[j], 0);
          }
        }
      }
      if(k==2){
        if(d==2){
          out[0] -= 0.5*(out[1]+out[3]);
          out[2] -= 0.5*(out[1]+out[4]);
          out[5] -= 0.5*(out[3]+out[6]);
          out[7] -= 0.5*(out[4]+out[6]);
        }
        if(d==3){
          out[0]  -= 0.5*(out[1] +out[3] +out[8]);
          out[2]  -= 0.5*(out[1] +out[4] +out[9]);
          out[5]  -= 0.5*(out[3] +out[6] +out[10]);
          out[7]  -= 0.5*(out[4] +out[6] +out[11]);
          out[12] -= 0.5*(out[8] +out[13]+out[15]);
          out[14] -= 0.5*(out[9] +out[13]+out[16]);
          out[17] -= 0.5*(out[10]+out[15]+out[18]);
          out[19] -= 0.5*(out[11]+out[16]+out[18]);
        }
      }
 
      //========================Debug test===========================
      /*for(size_t h =0; h < 9; h++){
        Dune::FieldVector<double,d> inh = {(h%3)/2.0,((h/3)%3)/2.0};
        std::cout << "At point " << inh << std::endl;
 
        Dune::FieldVector<double,8> outh;
        DuneVector alpha(0);
        for (size_t i=0; i<size(); i++)
        {
          // convert index i to multiindex
          if(i>0){ multiindex(alpha); }
 
          outh[i] = 1.0;
          for(int j = 0; j < d; j++){
            if(isEdge(alpha))
              outh[i] *= edgeVal(alpha[j], inh[j], 0);
            else{ //isNode
              outh[i] *= nodeVal(alpha[j], inh[j], 0);
            }
          }
        }
        if(d==2){
          outh[0] = outh[0] - 0.5*(outh[1]+outh[3]);
          outh[2] = outh[2] - 0.5*(outh[1]+outh[4]);
          outh[5] = outh[5] - 0.5*(outh[3]+outh[6]);
          outh[7] = outh[7] - 0.5*(outh[4]+outh[6]);
        }
 
        for (size_t i=0; i<size(); i++)
          std::cout << " value " << i << " is " << outh[i] << std::endl;
      }*/
      //===================End======================================
 
    }
 
    inline void
      evaluateJacobian (const typename Traits::DomainType& in,
          std::vector<typename Traits::JacobianType>& out) const
      {
        out.resize(size());
        DuneVector alpha(0);
 
        // Loop over all shape functions
        for (size_t i=0; i<size(); i++)
        {
          // convert index i to multiindex
          if(i>0){ multiindex(alpha); }
 
          out[i] =1.0;
          // 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
            if(isEdge(alpha))
              out[i][0][j] *= edgeVal(alpha[j],in[j], 1);
            else
              out[i][0][j] *= nodeVal(alpha[j],in[j], 1);
            // rest of the product
            for (int l=0; l<d; l++){
              if (l!=j){
                if(isEdge(alpha))
                  out[i][0][j] *= edgeVal(alpha[l],in[l],0);
                else
                  out[i][0][j] *= nodeVal(alpha[l],in[l],0);
              }
            }
          }
        }
        if(k==2){
        if(d==2){
          for(int j = 0; j < d; j++){
            out[0][0][j] -= 0.5*(out[1][0][j]+out[3][0][j]);
            out[2][0][j] -= 0.5*(out[1][0][j]+out[4][0][j]);
            out[5][0][j] -= 0.5*(out[3][0][j]+out[6][0][j]);
            out[7][0][j] -= 0.5*(out[4][0][j]+out[6][0][j]);
          }
        }
        if(d==3){
          for(int j = 0; j < d; j++){
            out[0][0][j]  -= 0.5*(out[1][0][j] +out[3][0][j] +out[8][0][j]);
            out[2][0][j]  -= 0.5*(out[1][0][j] +out[4][0][j] +out[9][0][j]);
            out[5][0][j]  -= 0.5*(out[3][0][j] +out[6][0][j] +out[10][0][j]);
            out[7][0][j]  -= 0.5*(out[4][0][j] +out[6][0][j] +out[11][0][j]);
            out[12][0][j] -= 0.5*(out[8][0][j] +out[13][0][j]+out[15][0][j]);
            out[14][0][j] -= 0.5*(out[9][0][j] +out[13][0][j]+out[16][0][j]);
            out[17][0][j] -= 0.5*(out[10][0][j]+out[15][0][j]+out[18][0][j]);
            out[19][0][j] -= 0.5*(out[11][0][j]+out[16][0][j]+out[18][0][j]);
          }
        }
        }
      }
 
    template<int diffOrder>
    inline void evaluate(
      const std::array<int,1>& direction,
      const typename Traits::DomainType& in,
      std::vector<typename Traits::RangeType>& out) const
    {
      static_assert(diffOrder == 1, "We only can compute first derivatives");
      out.resize(size());
      DuneVector alpha(0);
 
      // Loop over all shape functions
      for (size_t i=0; i<size(); i++)
      {
        // convert index i to multiindex
        if(i>0){ multiindex(alpha);}
 
        // Loop over all coordinate directions
        std::size_t j = direction[0];
 
        out[i] =1.0;
        // Initialize: the overall expression is a product
        // if j-th bit of i is set to -1, else 1
          if(isEdge(alpha))
            out[i][0] *= edgeVal(alpha[j],in[j],1);
          else
            out[i][0] *= nodeVal(alpha[j],in[j],1);
 
        // rest of the product
        for (std::size_t l=0; l<d; l++){
          if (l!=j){
            if(isEdge(alpha))
              out[i][0] *= edgeVal(alpha[l],in[l],0);
            else
              out[i][0] *= nodeVal(alpha[l],in[l],0);
          }
        }
      }
      if(k==2){
      if(d==2){
        out[0][0] -= 0.5*(out[1][0]+out[3][0]);
        out[2][0] -= 0.5*(out[1][0]+out[4][0]);
        out[5][0] -= 0.5*(out[3][0]+out[6][0]);
        out[7][0] -= 0.5*(out[4][0]+out[6][0]);
      }
      if(d==3){
        out[0][0]  -= 0.5*(out[1][0] +out[3][0] +out[8][0]);
        out[2][0]  -= 0.5*(out[1][0] +out[4][0] +out[9][0]);
        out[5][0]  -= 0.5*(out[3][0] +out[6][0] +out[10][0]);
        out[7][0]  -= 0.5*(out[4][0] +out[6][0] +out[11][0]);
        out[12][0] -= 0.5*(out[8][0] +out[13][0]+out[15][0]);
        out[14][0] -= 0.5*(out[9][0] +out[13][0]+out[16][0]);
        out[17][0] -= 0.5*(out[10][0]+out[15][0]+out[18][0]);
        out[19][0] -= 0.5*(out[11][0]+out[16][0]+out[18][0]);
      }
      }
    }
 
    unsigned int order () const
    {
      return k;
    }
  };
}
 
#endif