monomiallocalbasis.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_LOCALFUNCTIONS_MONOMIAL_MONOMIALLOCALBASIS_HH
#define DUNE_LOCALFUNCTIONS_MONOMIAL_MONOMIALLOCALBASIS_HH
 
#include <array>
#include <cassert>
#include <numeric>
 
#include <dune/common/fmatrix.hh>
 
#include "../common/localbasis.hh"
 
namespace Dune
{
  namespace MonomImp {
    template<int d, int k>
    struct Size {
      enum { val = Size<d,k-1>::val+Size<d-1,k>::val };
    };
    template<int d>
    struct Size<d, 0> {
      enum { val = 1 };
    };
    template<int k>
    struct Size<0, k> {
      enum { val = 1 };
    };
    template<>
    struct Size<0, 0> {
      enum { val = 1 };
    };
 
    template<class T>
    T ipow(T base, int exp)
    {
      T result(1);
      while (exp)
      {
        if (exp & 1)
          result *= base;
        exp >>= 1;
        base *= base;
      }
      return result;
    }
 
    template <typename Traits>
    class EvalAccess {
      std::vector<typename Traits::RangeType> &out;
#ifndef NDEBUG
      unsigned int first_unused_index;
#endif
 
    public:
      EvalAccess(std::vector<typename Traits::RangeType> &out_)
        : out(out_)
#ifndef NDEBUG
          , first_unused_index(0)
#endif
      { }
#ifndef NDEBUG
      ~EvalAccess() {
        assert(first_unused_index == out.size());
      }
#endif
      typename Traits::RangeFieldType &operator[](unsigned int index)
      {
        assert(index < out.size());
#ifndef NDEBUG
        if(first_unused_index <= index)
          first_unused_index = index+1;
#endif
        return out[index][0];
      }
    };
 
    template <typename Traits>
    class JacobianAccess {
      std::vector<typename Traits::JacobianType> &out;
      unsigned int row;
#ifndef NDEBUG
      unsigned int first_unused_index;
#endif
 
    public:
      JacobianAccess(std::vector<typename Traits::JacobianType> &out_,
                     unsigned int row_)
        : out(out_), row(row_)
#ifndef NDEBUG
          , first_unused_index(0)
#endif
      { }
#ifndef NDEBUG
      ~JacobianAccess() {
        assert(first_unused_index == out.size());
      }
#endif
      typename Traits::RangeFieldType &operator[](unsigned int index)
      {
        assert(index < out.size());
#ifndef NDEBUG
        if(first_unused_index <= index)
          first_unused_index = index+1;
#endif
        return out[index][0][row];
      }
    };
 
    template <typename Traits, int c>
    struct Evaluate
    {
      enum {
        d = Traits::dimDomain - c
      };
      template <typename Access>
      static void eval ( 
        const typename Traits::DomainType &in,
        const std::array<int, Traits::dimDomain> &derivatives,
        typename Traits::RangeFieldType prod,
        int bound,
        int& index,
        Access &access)
      {
        // start with the highest exponent for this dimension, then work down
        for (int e = bound; e >= 0; --e)
        {
          // the rest of the available exponents, to be used by the other
          // dimensions
          int newbound = bound - e;
          if(e < derivatives[d])
            Evaluate<Traits,c-1>::
            eval(in, derivatives, 0, newbound, index, access);
          else {
            int coeff = 1;
            for(int i = e - derivatives[d] + 1; i <= e; ++i)
              coeff *= i;
            // call the evaluator for the next dimension
            Evaluate<Traits,c-1>::
            eval(  // pass the coordinate and the derivatives unchanged
              in, derivatives,
              // also pass the product accumulated so far, but also
              // include the current dimension
              prod * ipow(in[d], e-derivatives[d]) * coeff,
              // pass the number of remaining exponents to the next
              // dimension
              newbound,
              // pass the next index to fill and the output access
              // wrapper
              index, access);
          }
        }
      }
    };
 
    template <typename Traits>
    struct Evaluate<Traits, 1>
    {
      enum { d = Traits::dimDomain-1 };
      template <typename Access>
      static void eval (const typename Traits::DomainType &in,
                        const std::array<int, Traits::dimDomain> &derivatives,
                        typename Traits::RangeFieldType prod,
                        int bound, int& index, Access &access)
      {
        if(bound < derivatives[d])
          prod = 0;
        else {
          int coeff = 1;
          for(int i = bound - derivatives[d] + 1; i <= bound; ++i)
            coeff *= i;
          prod *= ipow(in[d], bound-derivatives[d]) * coeff;
        }
        access[index] = prod;
        ++index;
      }
    };
 
  } //namespace MonomImp
 
  template<class D, class R, unsigned int d, unsigned int p>
  class MonomialLocalBasis
  {
  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 MonomImp::Size<d,p>::val;
    }
 
    inline void evaluateFunction (const typename Traits::DomainType& in,
                                  std::vector<typename Traits::RangeType>& out) const
    {
      out.resize(size());
      int index = 0;
      std::array<int, d> derivatives;
      std::fill(derivatives.begin(), derivatives.end(), 0);
      MonomImp::EvalAccess<Traits> access(out);
      for (unsigned int lp = 0; lp <= p; ++lp)
        MonomImp::Evaluate<Traits, d>::eval(in, derivatives, 1, lp, index, access);
    }
 
    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());
      int index = 0;
      std::array<int, d> derivatives;  // We need 'order' array as a signed int array
      std::copy(order.begin(), order.end(), derivatives.begin());
      MonomImp::EvalAccess<Traits> access(out);
      for (unsigned int lp = 0; lp <= p; ++lp)
        MonomImp::Evaluate<Traits, d>::eval(in, derivatives, 1, lp, index, access);
    }
 
    inline void
    evaluateJacobian (const typename Traits::DomainType& in,         // position
                      std::vector<typename Traits::JacobianType>& out) const      // return value
    {
      out.resize(size());
      std::array<int, d> derivatives;
      for(unsigned int i = 0; i < d; ++i)
        derivatives[i] = 0;
      for(unsigned int i = 0; i < d; ++i)
      {
        derivatives[i] = 1;
        int index = 0;
        MonomImp::JacobianAccess<Traits> access(out, i);
        for(unsigned int lp = 0; lp <= p; ++lp)
          MonomImp::Evaluate<Traits, d>::eval(in, derivatives, 1, lp, index, access);
        derivatives[i] = 0;
      }
    }
 
    unsigned int order () const
    {
      return p;
    }
  };
 
}
 
#endif // DUNE_LOCALFUNCTIONS_MONOMIAL_MONOMIALLOCALBASIS_HH