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_WHITNEY_EDGES0_5_INTERPOLATION_HH
#define DUNE_LOCALFUNCTIONS_WHITNEY_EDGES0_5_INTERPOLATION_HH
 
#include <cstddef>
#include <vector>
 
#include <dune/localfunctions/whitney/edges0.5/common.hh>
 
namespace Dune {
 
  //
  // Interpolation
  //
 
 
  template<class Geometry, class Traits_>
  class EdgeS0_5Interpolation :
    private EdgeS0_5Common<Traits_::dimDomainLocal,
        typename Traits_::DomainField>
  {
  public:
    typedef Traits_ Traits;
 
  private:
    static const std::size_t dim = Traits::dimDomainLocal;
    typedef EdgeS0_5Common<dim, typename Traits::DomainField> Base;
    using Base::refelem;
    using Base::s;
 
    std::vector<typename Traits::DomainGlobal> edgev;
 
  public:
 
    template<typename VertexOrder>
    EdgeS0_5Interpolation(const Geometry& geo,
                          const VertexOrder& vertexOrder) :
      edgev(s)
    {
      for(std::size_t i = 0; i < s; ++i) {
        const std::size_t i0 = refelem.subEntity(i,dim-1,0,dim);
        const std::size_t i1 = refelem.subEntity(i,dim-1,1,dim);
 
        edgev[i] = geo.corner(i1);
        edgev[i] -= geo.corner(i0);
        edgev[i] /= edgev[i].two_norm();
 
        const typename VertexOrder::iterator& edgeVertexOrder =
          vertexOrder.begin(dim-1, i);
        if(edgeVertexOrder[0] > edgeVertexOrder[1])
          edgev[i] *= -1;
      }
    }
 
    template<typename F, typename C>
    void interpolate(const F& f, std::vector<C>& out) const {
      typename Traits::Range y;
 
      out.resize(s);
 
      for(std::size_t i = 0; i < s; ++i) {
        f.evaluate(refelem.position(i,dim-1), y);
 
        out[i] = y * edgev[i];
      }
    }
  };
 
} // namespace Dune
 
#endif // DUNE_LOCALFUNCTIONS_WHITNEY_EDGES0_5_INTERPOLATION_HH