interpolation.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_INTERPOLATION_HH
#define DUNE_LOCALFUNCTIONS_META_POWER_INTERPOLATION_HH
 
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <vector>
#include <dune/common/function.hh>
#include <dune/localfunctions/common/localinterpolation.hh>
 
namespace Dune {
 
 
  template<class Backend, class BasisTraits>
  class PowerInterpolation {
    static_assert(Backend::Traits::dimRange == 1,
                  "PowerInterpolation  works only with scalar backends");
 
    const Backend *backend;
 
  public:
    typedef BasisTraits Traits;
 
 
    PowerInterpolation(const Backend &backend_) : backend(&backend_) { }
 
  private:
    template<class F>
    class ComponentEvaluator :
      public Function<typename Backend::Traits::DomainLocal, typename Backend::Traits::Range>
    {
      const F &f;
      std::size_t comp;
 
    public:
      ComponentEvaluator(const F &f_, std::size_t comp_) :
        f(f_), comp(comp_)
      { }
 
      void evaluate(const typename Backend::Traits::DomainLocal &x,
                    typename Backend::Traits::Range &y) const
      {
        typename Traits::Range fy = f(x);
        y[0] = fy[comp];
      }
    };
 
  public:
 
    template<typename F, typename C>
    void interpolate(const F& ff, std::vector<C>& out) const {
 
      auto&& f = Impl::makeFunctionWithCallOperator<typename Backend::Traits::DomainLocal>(ff);
 
 
      out.clear();
      std::vector<C> cout;
      for(std::size_t d = 0; d < Traits::dimRange; ++d) {
        // When dropping support for `evaluate()` we can simply use a lambda
        // instead of ComponentEvaluator. But changing this now would break
        // PowerInterpolation for FE-implementation outside of dune-localfunctions
        // which may not have been adjusted so far.
        backend->interpolate(ComponentEvaluator<std::decay_t<decltype(f)>>(f, d), cout);
        if(d == 0)
          out.resize(cout.size()*Traits::dimRange);
        // make sure the size of cout does not change surprisingly
        assert(out.size() == cout.size()*Traits::dimRange);
        std::copy(cout.begin(), cout.end(), out.begin() + d*cout.size());
      }
    }
  };
 
} // namespace Dune
 
#endif // DUNE_LOCALFUNCTIONS_META_POWER_INTERPOLATION_HH