qkdglobatto.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_PDELAB_FINITEELEMENT_QKDGLOBATTO_HH
#define DUNE_PDELAB_FINITEELEMENT_QKDGLOBATTO_HH
 
#include <numeric>
 
#include <dune/localfunctions/common/localbasis.hh>
#include <dune/localfunctions/common/localkey.hh>
#include <dune/localfunctions/common/localfiniteelementtraits.hh>
#include <dune/localfunctions/common/localtoglobaladaptors.hh>
#include <dune/pdelab/finiteelement/qkdglagrange.hh>
 
namespace Dune
{
  namespace QkStuff
  {
 
    template<class D, class R, int k>
    class GaussLobattoLagrangePolynomials
    {
      R xi_gl[k+1];
      R w_gl[k+1];
 
    public:
      GaussLobattoLagrangePolynomials ()
      {
        int matched_order=-1;
        int matched_size=-1;
        for (int order=1; order<=40; order++)
          {
            const Dune::QuadratureRule<D,1>& rule = Dune::QuadratureRules<D,1>::rule(Dune::GeometryType::cube,order,Dune::QuadratureType::GaussLobatto);
            if (rule.size()==k+1)
              {
                matched_order = order;
                matched_size = rule.size();
                //std::cout << "GL: input order=" << order << " delivered=" << rule.order() << " size=" << rule.size()<< std::endl;
                break;
              }
          }
        if (matched_order<0) DUNE_THROW(Dune::Exception,"could not find Gauss Lobatto rule of appropriate size");
        const Dune::QuadratureRule<D,1>& rule = Dune::QuadratureRules<D,1>::rule(Dune::GeometryType::cube,matched_order,Dune::QuadratureType::GaussLobatto);
        size_t count=0;
        for (typename Dune::QuadratureRule<D,1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
          {
            size_t group=count/2;
            size_t member=count%2;
            size_t newj;
            if (member==1) newj=group; else newj=k-group;
            xi_gl[newj] = it->position()[0];
            w_gl[newj] = it->weight();
            count++;
          }
        for (int j=0; j<matched_size/2; j++)
          if (xi_gl[j]>0.5)
            {
              R temp=xi_gl[j];
              xi_gl[j] = xi_gl[k-j];
              xi_gl[k-j] = temp;
              temp=w_gl[j];
              w_gl[j] = w_gl[k-j];
              w_gl[k-j] = temp;
            }
        // for (int i=0; i<k+1; i++)
        //   std::cout << "i=" << i
        //             << " xi=" << xi_gl[i]
        //             << " w=" << w_gl[i]
        //             << std::endl;
      }
 
      // ith Lagrange polynomial of degree k in one dimension
      R p (int i, D x) const
      {
        R result(1.0);
        for (int j=0; j<=k; j++)
          if (j!=i) result *= (x-xi_gl[j])/(xi_gl[i]-xi_gl[j]);
        return result;
      }
 
      // derivative of ith Lagrange polynomial of degree k in one dimension
      R dp (int i, D x) const
      {
        R result(0.0);
 
        for (int j=0; j<=k; j++)
          if (j!=i)
            {
              R prod( 1.0/(xi_gl[i]-xi_gl[j]) );
              for (int l=0; l<=k; l++)
                if (l!=i && l!=j) prod *= (x-xi_gl[l])/(xi_gl[i]-xi_gl[l]);
              result += prod;
            }
        return result;
      }
 
      // get ith Lagrange point
      R x (int i) const
      {
        return xi_gl[i];
      }
 
      // get weight of ith Lagrange point
      R w (int i) const
      {
        return w_gl[i];
      }
    };
 
    template<class D, class R, int k, int d>
    class QkGLLocalBasis
    {
      enum{ n = QkSize<k,d>::value };
      GaussLobattoLagrangePolynomials<D,R,k> poly;
 
    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 QkSize<k,d>::value;
      }
 
      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<k,d>(i));
 
            // initialize product
            out[i] = 1.0;
 
            // dimension by dimension
            for (int j=0; j<d; j++)
              out[i] *= poly.p(alpha[j],in[j]);
          }
      }
 
      inline void
      evaluateJacobian (const typename Traits::DomainType& in,         // position
                        std::vector<typename Traits::JacobianType>& out) const      // return value
      {
        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<k,d>(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] = poly.dp(alpha[j],in[j]);
 
                // rest of the product
                for (int l=0; l<d; l++)
                  if (l!=j)
                    out[i][0][j] *= poly.p(alpha[l],in[l]);
              }
          }
      }
 
      void partial(const std::array<unsigned int, Traits::dimDomain>& order,
                   const typename Traits::DomainType& in,
                   std::vector<typename Traits::RangeType>& out) const {
        auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
        if (totalOrder == 0) {
          evaluateFunction(in, out);
        } else {
          DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
        }
      }
 
      unsigned int order () const
      {
        return k;
      }
    };
 
    template<int k, int d, class LB>
    class QkGLLocalInterpolation
    {
      GaussLobattoLagrangePolynomials<double,double,k> poly;
 
    public:
 
      template<typename F, typename C>
      void interpolate (const F& f, std::vector<C>& out) const
      {
        typename LB::Traits::DomainType x;
 
        out.resize(QkSize<k,d>::value);
 
        for (int i=0; i<QkSize<k,d>::value; i++)
          {
            // convert index i to multiindex
            Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
 
            // Generate coordinate of the i-th Lagrange point
            for (int j=0; j<d; j++)
              x[j] = poly.x(alpha[j]);
 
            out[i] = f(x);
          }
      }
    };
 
    template<int d, class LB>
    class QkGLLocalInterpolation<0,d,LB>
    {
    public:
      template<typename F, typename C>
      void interpolate (const F& f, std::vector<C>& out) const
      {
        typename LB::Traits::DomainType x(0.5);
        out.resize(1);
        out[0] = f(x);
      }
    };
 
  }
 
  template<class D, class R, int k, int d>
  class QkDGGLLocalFiniteElement
  {
    typedef QkStuff::QkGLLocalBasis<D,R,k,d> LocalBasis;
    typedef QkStuff::QkDGLocalCoefficients<k,d> LocalCoefficients;
    typedef QkStuff::QkGLLocalInterpolation<k,d,LocalBasis> LocalInterpolation;
 
  public:
    // static number of basis functions
    enum{ n = QkStuff::QkSize<k,d>::value };
 
    typedef LocalFiniteElementTraits<LocalBasis,LocalCoefficients,LocalInterpolation> Traits;
 
    QkDGGLLocalFiniteElement ()
      : gt(GeometryTypes::cube(d))
    {}
 
    const typename Traits::LocalBasisType& localBasis () const
    {
      return basis;
    }
 
    const typename Traits::LocalCoefficientsType& localCoefficients () const
    {
      return coefficients;
    }
 
    const typename Traits::LocalInterpolationType& localInterpolation () const
    {
      return interpolation;
    }
 
    std::size_t size() const
    {
      return basis.size();
    }
 
    GeometryType type () const
    {
      return gt;
    }
 
    QkDGGLLocalFiniteElement* clone () const
    {
      return new QkDGGLLocalFiniteElement(*this);
    }
 
  private:
    LocalBasis basis;
    LocalCoefficients coefficients;
    LocalInterpolation interpolation;
    GeometryType gt;
  };
 
 
  template<class Geometry, class RF, int k>
  class QkDGGLFiniteElementFactory :
    public ScalarLocalToGlobalFiniteElementAdaptorFactory<
    QkDGGLLocalFiniteElement<
      typename Geometry::ctype, RF, k, Geometry::mydimension
      >,
    Geometry
    >
  {
    typedef QkDGGLLocalFiniteElement<
      typename Geometry::ctype, RF, k, Geometry::mydimension
      > LFE;
    typedef ScalarLocalToGlobalFiniteElementAdaptorFactory<LFE, Geometry> Base;
 
    static const LFE lfe;
 
  public:
    QkDGGLFiniteElementFactory() : Base(lfe) {}
  };
 
  template<class Geometry, class RF, int k>
  const typename QkDGGLFiniteElementFactory<Geometry, RF, k>::LFE
  QkDGGLFiniteElementFactory<Geometry, RF, k>::lfe;
 
}
 
#endif // DUNE_PDELAB_FINITEELEMENT_QKDGLOBATTO_HH