basis.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_LOCALFUNCTIONS_META_POWER_BASIS_HH
#define DUNE_LOCALFUNCTIONS_META_POWER_BASIS_HH
 
#include <numeric>
#include <cstddef>
#include <vector>
 
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
 
namespace Dune {
 
 
  template<class Backend, std::size_t dimR>
  class PowerBasis {
    static_assert(Backend::Traits::dimRange == 1,
                  "PowerBasis works only with scalar backends");
 
    // don't use a reference here so this class stays copyable
    const Backend *backend;
 
  public:
    struct Traits : public Backend::Traits
    {
      static const std::size_t dimRange = dimR;
      typedef FieldVector<typename Traits::RangeField, dimR> Range;
 
 
      typedef FieldMatrix<typename Traits::RangeField, dimR,
          Traits::dimDomainGlobal> Jacobian;
    };
 
 
    PowerBasis(const Backend &backend_) : backend(&backend_) { }
 
    std::size_t size () const { return backend->size()*dimR; }
    std::size_t order () const { return backend->order(); }
 
    void evaluateFunction(const typename Traits::DomainLocal& in,
                          std::vector<typename Traits::Range>& out) const
    {
      std::vector<typename Backend::Traits::Range> backendValues;
      backend->evaluateFunction(in, backendValues);
      out.assign(size(), typename Traits::Range(0));
      for(std::size_t d = 0; d < dimR; ++d)
        for(std::size_t i = 0; i < backend->size(); ++i)
          out[d*backend->size()+i][d] = backendValues[i][0];
    }
 
    void evaluateJacobian(const typename Traits::DomainLocal& in,
                          std::vector<typename Traits::Jacobian>& out) const
    {
      std::vector<typename Backend::Traits::Jacobian> backendValues;
      backend->evaluateJacobian(in, backendValues);
      out.assign(size(), typename Traits::Jacobian(0));
      for(std::size_t d = 0; d < dimR; ++d)
        for(std::size_t i = 0; i < backend->size(); ++i)
          out[d*backend->size()+i][d] = backendValues[i][0];
    }
 
    void partial (const std::array<unsigned int, Backend::Traits::dimDomainGlobal>& order,
                  const typename Traits::DomainLocal& in,         // position
                  std::vector<typename Traits::Range>& out) const      // return value
    {
      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");
      }
    }
  };
 
} // namespace Dune
 
#endif // DUNE_LOCALFUNCTIONS_META_POWER_BASIS_HH