function.hh

//-*- tab-width: 4; indent-tabs-mode: nil -*-
#ifndef DUNE_PDELAB_FUNCTION_HH
#define DUNE_PDELAB_FUNCTION_HH
 
#include <iostream>
#include <sstream>
 
#include <dune/common/deprecated.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
 
#include <dune/grid/utility/hierarchicsearch.hh>
 
#include <dune/typetree/typetree.hh>
 
#include "vtkexport.hh"
#include "geometrywrapper.hh"
#include "typetraits.hh"
 
namespace Dune {
  namespace PDELab {
 
 
    template<class DF, int n, class D, class RF, int m, class R>
    struct FunctionTraits
    {
      typedef DF DomainFieldType;
 
      enum {
        dimDomain = n
      };
 
      typedef D DomainType;
 
      typedef RF RangeFieldType;
 
      enum {
        dimRange = m
      };
 
      typedef R RangeType;
    };
 
    template<class T, class Imp>
    class FunctionInterface
    {
    public:
      typedef T Traits;
 
      inline void evaluate (const typename Traits::DomainType& x,
                            typename Traits::RangeType& y) const
      {
        asImp().evaluate(x,y);
      }
 
    private:
      Imp& asImp () {return static_cast<Imp &> (*this);}
      const Imp& asImp () const {return static_cast<const Imp &>(*this);}
    };
 
    class InstationaryFunctionDefaults
    {
    public:
 
      template<typename Time>
      inline void setTime(Time t)
      { }
    };
 
    template<typename GV>
    struct PowerCompositeGridFunctionTraits
    {
      typedef GV GridViewType;
 
      typedef typename GV::Traits::template Codim<0>::Entity ElementType;
 
    };
 
 
    class GridFunctionOutputParameters
    {
 
    public:
 
      struct Output
      {
 
        enum DataSetType
          {
            vertexData, 
            cellData    
          };
      };
 
 
      GridFunctionOutputParameters(Output::DataSetType dataSetType = Output::vertexData)
        : _dataSetType(dataSetType)
      {}
 
      Output::DataSetType dataSetType() const
      {
        return _dataSetType;
      }
 
      void setDataSetType(Output::DataSetType dataSetType)
      {
        _dataSetType = dataSetType;
      }
 
    private:
 
      Output::DataSetType _dataSetType;
 
    };
 
    template<class GV, class RF, int m, class R>
    struct GridFunctionTraits
      : public FunctionTraits<typename GV::Grid::ctype, GV::dimension,
                              Dune::FieldVector<typename GV::Grid::ctype,
                                                GV::dimension>,
                              RF, m, R>
      , public PowerCompositeGridFunctionTraits<GV>
    {
    };
 
    template<class T, class Imp>
    class GridFunctionInterface
      : public GridFunctionOutputParameters
    {
    public:
      typedef T Traits;
 
      GridFunctionInterface(Output::DataSetType dataSetType = Output::vertexData)
        : GridFunctionOutputParameters(dataSetType)
      {}
 
      inline void evaluate (const typename Traits::ElementType& e,
                            const typename Traits::DomainType& x,
                            typename Traits::RangeType& y) const
      {
        asImp().evaluate(e,x,y);
      }
 
      inline const typename Traits::GridViewType& getGridView () const
      {
        return asImp().getGridView();
      }
 
    private:
      Imp& asImp () {return static_cast<Imp &> (*this);}
      const Imp& asImp () const {return static_cast<const Imp &>(*this);}
    };
 
    template<class GV, class RF, int m, class R>
    struct BoundaryGridFunctionTraits
      : public FunctionTraits<typename GV::Grid::ctype, GV::dimension-1,
                              Dune::FieldVector<typename GV::Grid::ctype,
                                                GV::dimension-1>,
                              RF, m, R>
    {
      typedef GV GridViewType;
    };
 
 
    // \tparam T The type of the BoundaryGridFunctionTraits.
    // \tparam Imp The type of the implementing class.
    template<class T, class Imp>
    class BoundaryGridFunctionInterface
    {
    public:
      typedef T Traits;
 
      template<typename I>
      inline void evaluate (const IntersectionGeometry<I>& ig,
                            const typename Traits::DomainType& x,
                            typename Traits::RangeType& y) const
      {
        asImp().evaluate(ig,x,y);
      }
 
      inline const typename Traits::GridViewType& getGridView () const
      {
        return asImp().getGridView();
      }
 
    private:
      Imp& asImp () {return static_cast<Imp &> (*this);}
      const Imp& asImp () const {return static_cast<const Imp &>(*this);}
    };
 
    //============================
    // Function tree
    //============================
 
 
    struct GridFunctionTag {};
 
    template<class T, class Imp>
    class GridFunctionBase
      : public GridFunctionInterface<T,Imp>
      , public TypeTree::LeafNode
    {
      using Base = GridFunctionInterface<T,Imp>;
    public:
      typedef GridFunctionTag ImplementationTag;
      typedef typename T::GridViewType GridViewType;
 
      using Output = typename Base::Output;
 
      GridFunctionBase(typename Output::DataSetType dataSetType = Output::vertexData)
        : Base(dataSetType)
      {}
    };
 
 
    template<class T, class Imp>
    class BoundaryGridFunctionBase
      : public BoundaryGridFunctionInterface<T,Imp>
      , public TypeTree::LeafNode
    {
    public:
      typedef GridFunctionTag ImplementationTag;
      typedef typename T::GridViewType GridViewType;
    };
 
 
    template<typename TT>
    struct PowerCompositeSetTimeVisitor
      : public TypeTree::TreeVisitor, public TypeTree::DynamicTraversal
    {
      TT time;
      PowerCompositeSetTimeVisitor(const TT time_) : time(time_) {}
 
      template<typename LeafNode, typename TreePath>
      void leaf(LeafNode& node, TreePath treePath) const
      {
        node.setTime(time);
      }
    };
 
    struct PowerGridFunctionTag {};
 
    template<class T, std::size_t k>
    class PowerGridFunction
      : public TypeTree::PowerNode<T,k>
    {
 
      typedef TypeTree::PowerNode<T,k> BaseT;
 
    public:
 
      typedef PowerCompositeGridFunctionTraits<typename T::GridViewType> Traits;
 
      typedef PowerGridFunctionTag ImplementationTag;
 
      typedef typename T::GridViewType GridViewType;
 
      template <typename TT>
      void setTime(TT time){
        PowerCompositeSetTimeVisitor<TT> visitor(time);
        TypeTree::applyToTree(*this,visitor);
      }
 
      PowerGridFunction()
      {}
 
      PowerGridFunction (const std::array<std::shared_ptr<T>,k>& container)
        : BaseT(container) {}
 
      PowerGridFunction (T& t)
        : BaseT(t) {}
 
      // TODO: PowerGridFunction (T** t) : ...
 
#ifdef DOXYGEN
      PowerGridFunction (T& t0, T& t1, ...)
      {
      }
 
#else
 
      PowerGridFunction (T& c0,
                         T& c1)
        : BaseT(c0,c1)
      {
      }
 
      PowerGridFunction (T& c0,
                         T& c1,
                         T& c2)
        : BaseT(c0,c1,c2)
      {
      }
 
      PowerGridFunction (T& c0,
                         T& c1,
                         T& c2,
                         T& c3)
        : BaseT(c0,c1,c2,c3)
      {
      }
 
      PowerGridFunction (T& c0,
                         T& c1,
                         T& c2,
                         T& c3,
                         T& c4)
        : BaseT(c0,c1,c2,c3,c4)
      {
      }
 
      PowerGridFunction (T& c0,
                         T& c1,
                         T& c2,
                         T& c3,
                         T& c4,
                         T& c5)
        : BaseT(c0,c1,c2,c3,c4,c5)
      {
      }
 
      PowerGridFunction (T& c0,
                         T& c1,
                         T& c2,
                         T& c3,
                         T& c4,
                         T& c5,
                         T& c6)
        : BaseT(c0,c1,c2,c3,c4,c5,c6)
      {
      }
 
      PowerGridFunction (T& c0,
                         T& c1,
                         T& c2,
                         T& c3,
                         T& c4,
                         T& c5,
                         T& c6,
                         T& c7)
        : BaseT(c0,c1,c2,c3,c4,c5,c6,c7)
      {
      }
 
      PowerGridFunction (T& c0,
                         T& c1,
                         T& c2,
                         T& c3,
                         T& c4,
                         T& c5,
                         T& c6,
                         T& c7,
                         T& c8)
        : BaseT(c0,c1,c2,c3,c4,c5,c6,c7,c8)
      {
      }
 
      PowerGridFunction (T& c0,
                         T& c1,
                         T& c2,
                         T& c3,
                         T& c4,
                         T& c5,
                         T& c6,
                         T& c7,
                         T& c8,
                         T& c9)
        : BaseT(c0,c1,c2,c3,c4,c5,c6,c7,c8,c9)
      {
      }
 
#endif // DOXYGEN
    };
 
    struct CompositeGridFunctionTag {};
 
    template<typename... Children>
    class CompositeGridFunction
      : public TypeTree::CompositeNode<Children...>
    {
 
      typedef TypeTree::CompositeNode<Children...> BaseT;
 
    public:
 
      typedef CompositeGridFunctionTag ImplementationTag;
 
      typedef PowerCompositeGridFunctionTraits<typename BaseT::template Child<0>::Type::GridViewType> Traits;
 
      typedef typename BaseT::template Child<0>::Type::GridViewType GridViewType;
 
      CompositeGridFunction()
      {}
 
      CompositeGridFunction (Children&... children)
        : BaseT(TypeTree::assertGridViewType<typename BaseT::template Child<0>::Type>(children)...)
      {
      }
 
      CompositeGridFunction (std::shared_ptr<Children>... children)
        : BaseT(children...)
      {}
 
      template <typename TT>
      void setTime(TT time){
        PowerCompositeSetTimeVisitor<TT> visitor(time);
        TypeTree::applyToTree(*this,visitor);
      }
 
#ifdef DOXYGEN
      CompositeGridFunction (T0& t0, T1& t1, ...) {}
#endif //DOXYGEN
    };
 
    //========================================================
    // helper template to turn an ordinary GridFunction into a
    // GridFunctionTree leaf
    //========================================================
 
    template<class Imp>
    class GridFunctionBaseAdapter
      : public GridFunctionBase<typename Imp::Traits,
                                GridFunctionBaseAdapter<Imp> >
    {
      const Imp &imp;
 
    public:
 
      GridFunctionBaseAdapter(const Imp& imp_)
        : imp(imp_)
      { }
 
 
      inline void evaluate (const typename Imp::Traits::ElementType& e,
                            const typename Imp::Traits::DomainType& x,
                            typename Imp::Traits::RangeType& y) const
      {
        imp.evaluate(e,x,y);
      }
 
      inline const typename Imp::Traits::GridViewType& getGridView () const
      {
        return imp.getGridView();
      }
    };
 
    //=======================================
    // helper template for analytic functions
    //=======================================
 
    template<typename GV, typename RF, int m>
    struct AnalyticGridFunctionTraits
      : public GridFunctionTraits<GV, RF, m, Dune::FieldVector<RF,m> >
    {
    };
 
    template<typename T, typename Imp>
    class AnalyticGridFunctionBase
      : public GridFunctionBase<T,AnalyticGridFunctionBase<T,Imp> >
    {
    public:
      typedef T Traits;
 
      AnalyticGridFunctionBase (const typename Traits::GridViewType& g_) : g(g_) {}
 
      inline void evaluate (const typename Traits::ElementType& e,
                            const typename Traits::DomainType& x,
                            typename Traits::RangeType& y) const
      {
        asImp().evaluateGlobal(e.geometry().global(x),y);
      }
 
      inline const typename Traits::GridViewType& getGridView () const
      {
        return g;
      }
 
    private:
      typename Traits::GridViewType g;
      Imp& asImp () {return static_cast<Imp &> (*this);}
      const Imp& asImp () const {return static_cast<const Imp &>(*this);}
    };
 
 
    // Adapter takes a vector-valued grid function and provides evaluation
    // of normal flux on the interior of faces.
    template<typename T>
    class NormalFluxGridFunctionAdapter
      : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<typename T::Traits::GridViewType,
                                                                                    typename T::Traits::RangeFieldType,
                                                                                    1,
                                                                                    Dune::FieldVector<typename T::Traits::RangeFieldType,1>
                                                                                    >,
                                                   NormalFluxGridFunctionAdapter<T> >
    {
    public:
      typedef Dune::PDELab::GridFunctionTraits<typename T::Traits::GridViewType,typename T::Traits::RangeFieldType,1,Dune::FieldVector<typename T::Traits::RangeFieldType,1> > Traits;
      typedef Dune::PDELab::GridFunctionInterface<Traits,NormalFluxGridFunctionAdapter<T> > BaseT;
 
      NormalFluxGridFunctionAdapter (const T& t_) : t(stackobject_to_shared_ptr(t_)) {}
 
 
      inline void evaluate (const typename Traits::ElementType& e,
                            const typename Traits::DomainType& x,
                            typename Traits::RangeType& y) const
      {
        // ensure correct size
        static_assert((static_cast<int>(T::Traits::GridViewType::dimension)==static_cast<int>(T::Traits::dimRange)),"number of components must equal dimension");
 
        // evaluate velocity
        typename T::Traits::RangeType v;
        t->evaluate(e,x,v);
 
        // implementation only handles triangles so far
        if (!e.geometry().type().isTriangle())
          DUNE_THROW(Dune::NotImplemented, "only implemented for triangles");
 
        // start and end corner in local numbering
        int n0, n1;
 
        typename Traits::DomainType nu;
 
        // determine outer unit normal
        if (std::abs(x[0])<1E-10)
          {
            // edge 1
            n0 = 2;
            n1 = 0;
 
            nu = e.geometry().corner(n1);
            nu -= e.geometry().corner(n0);
            typename Traits::DomainFieldType temp = nu[0];
            nu[0] = nu[1];
            nu[1] = -temp;
            nu /= nu.two_norm();
            y = v[0]*nu[0]+v[1]*nu[1];
            return;
          }
 
        if (std::abs(x[1])<1E-10)
          {
            // edge 2
            n0 = 0;
            n1 = 1;
 
            nu = e.geometry().corner(n1);
            nu -= e.geometry().corner(n0);
            typename Traits::DomainFieldType temp = nu[0];
            nu[0] = nu[1];
            nu[1] = -temp;
            nu /= nu.two_norm();
            y = v[0]*nu[0]+v[1]*nu[1];
            return;
          }
 
        if (std::abs(x[0]+x[1]-1.0)<1E-10)
          {
            // edge 0
            n0 = 1;
            n1 = 2;
 
            nu = e.geometry().corner(n1);
            nu -= e.geometry().corner(n0);
            typename Traits::DomainFieldType temp = nu[0];
            nu[0] = nu[1];
            nu[1] = -temp;
            nu /= nu.two_norm();
            y = v[0]*nu[0]+v[1]*nu[1];
            return;
          }
 
        DUNE_THROW(Dune::Exception, "x needs to be on an edge");
      }
 
      inline const typename Traits::GridViewType& getGridView () const
      {
        return t->getGridView();
      }
 
    private:
      std::shared_ptr<T const> t;
    };
 
    // Adapter takes a vector-valued grid function and applies
    // backward Piola transformation on each element
    template<typename T>
    class PiolaBackwardAdapter
      : public Dune::PDELab::GridFunctionBase<typename T::Traits,PiolaBackwardAdapter<T> >
    {
    public:
      typedef typename T::Traits::GridViewType GridViewType;
      typedef typename T::Traits Traits;
      typedef Dune::PDELab::GridFunctionBase<Traits,PiolaBackwardAdapter<T> > BaseT;
      // typedef GridFunctionTag ImplementationTag;
 
      PiolaBackwardAdapter (const T& t_) : t(stackobject_to_shared_ptr(t_)) {}
 
 
      inline void evaluate (const typename Traits::ElementType& e,
                            const typename Traits::DomainType& x,
                            typename Traits::RangeType& y) const
      {
        // evaluate velocity
        typename T::Traits::RangeType v;
        t->evaluate(e,x,v);
 
        // apply Piola transformation
        typename Traits::ElementType::Geometry::JacobianInverseTransposed
          J = e.geometry().jacobianInverseTransposed(x);
        y = 0;
        J.umtv(v,y);
        y *= e.geometry().integrationElement(x);
      }
 
      inline const typename Traits::GridViewType& getGridView () const
      {
        return t->getGridView();
      }
 
    private:
      std::shared_ptr<T const> t;
    };
 
 
    //==========================
    // template metaprograms
    //==========================
 
    namespace {
 
      template<typename VTKWriter>
      struct AddGridFunctionsToVTKWriter
        : public TypeTree::TreeVisitor
        , public TypeTree::DynamicTraversal
      {
 
        VTKWriter& w;
        const std::string s;
 
        AddGridFunctionsToVTKWriter(VTKWriter& w_, const std::string & s_) :
          w(w_), s(s_) {}
 
        template<typename T, typename TreePath>
        void leaf(const T& t, TreePath treePath) {
          std::stringstream name;
          name << s;
          for (std::size_t i=0; i < treePath.size(); ++i)
            name << "_" << treePath.element(i);
          w.addVertexData(std::make_shared< VTKGridFunctionAdapter<T> >(t,name.str()));
        }
      };
 
    } // anonymous namespace
 
    template<typename GV, typename T>
    void vtkwriter_tree_addvertexdata (Dune::VTKWriter<GV>& w, const T& t, std::string s = "data")
    {
      AddGridFunctionsToVTKWriter<Dune::VTKWriter<GV> > visitor(w,s);
      TypeTree::applyToTree(t,visitor);
    }
 
 
 
    template<typename G, typename T>
    class FunctionToGridFunctionAdapter :
      public GridFunctionBase<GridFunctionTraits<
                                     G,
                                     typename T::Traits::RangeFieldType,
                                     T::Traits::dimRange,
                                     typename T::Traits::RangeType>,
                                   FunctionToGridFunctionAdapter<G,T> >
    {
    public:
      typedef GridFunctionTraits<G,
                                 typename T::Traits::RangeFieldType,
                                 T::Traits::dimRange,
                                 typename T::Traits::RangeType> Traits;
      static_assert(
                    (std::is_same<typename T::Traits::DomainFieldType,
                     typename Traits::DomainFieldType>::value),
                    "GridView's and wrapped Functions DomainFieldType don't match");
      static_assert(
                    T::Traits::dimDomain==Traits::dimDomain,
                    "GridView's and wrapped Functions dimDomain don't match");
      static_assert(
                    (std::is_same<typename T::Traits::DomainType,
                     typename Traits::DomainType>::value),
                    "GridView's and wrapped Functions DomainType don't match");
 
      FunctionToGridFunctionAdapter (const G& g_, const T& t_) : g(g_), t(t_) {}
 
      inline void evaluate (const typename Traits::ElementType& e,
                            const typename Traits::DomainType& x,
                            typename Traits::RangeType& y) const
      {
        t.evaluate(e.geometry().global(x),y);
      }
 
      inline const typename Traits::GridViewType& getGridView () const
      {
        return g;
      }
 
    private:
      G g;
      const T& t;
    };
 
    template<typename GF>
    class GridFunctionToFunctionAdapter
      : public FunctionInterface<FunctionTraits<typename GF::Traits::GridViewType::ctype,
                                                GF::Traits::GridViewType::dimensionworld,
                                                Dune::FieldVector<typename GF::Traits::GridViewType::ctype,
                                                                  GF::Traits::GridViewType::dimensionworld
                                                                  >,
                                                typename GF::Traits::RangeFieldType,
                                                GF::Traits::dimRange,
                                                Dune::FieldVector<typename GF::Traits::RangeFieldType,
                                                                  GF::Traits::dimRange>
                                                >,
                                 GridFunctionToFunctionAdapter<GF> >
    {
    public:
      typedef FunctionTraits<typename GF::Traits::GridViewType::ctype,
                             GF::Traits::GridViewType::dimensionworld,
                             Dune::FieldVector<typename GF::Traits::GridViewType::ctype,
                                               GF::Traits::GridViewType::dimensionworld
                                               >,
                             typename GF::Traits::RangeFieldType,
                             GF::Traits::dimRange,
                             Dune::FieldVector<typename GF::Traits::RangeFieldType,
                                               GF::Traits::dimRange>
                             > Traits;
 
      GridFunctionToFunctionAdapter(const GF &gf_)
        : gf(gf_)
        , hsearch(gf.getGridView().grid(), gf.getGridView().indexSet())
      { }
 
      inline void evaluate (const typename Traits::DomainType& x,
                            typename Traits::RangeType& y) const
      {
        typename GF::Traits::GridViewType::Grid::Traits::template Codim<0>::Entity
          e = hsearch.findEntity(x);
        gf.evaluate(e, e.geometry().local(x), y);
      }
 
    private:
      const GF &gf;
      const Dune::HierarchicSearch<typename GF::Traits::GridViewType::Grid,
                                   typename GF::Traits::GridViewType::IndexSet> hsearch;
    };
 
 
    template<typename T, typename E>
    class GlobalFunctionToLocalFunctionAdapter :
      public FunctionInterface<typename T::Traits,
                               GlobalFunctionToLocalFunctionAdapter<T,E> >
    {
    public:
      typedef typename T::Traits Traits;
 
      GlobalFunctionToLocalFunctionAdapter (const T& t_, const E& e_) : t(t_), e(e_) {}
 
      inline void evaluate (const typename Traits::DomainType& x,
                            typename Traits::RangeType& y) const
      {
        t.evaluate(e.geometry().global(x),y);
      }
 
    private:
      const T& t;
      const E& e;
    };
 
 
    template<typename T> // T: GridFunction, E: Entity
    class GridFunctionToLocalFunctionAdapter :
      public FunctionInterface<typename T::Traits,
                               GridFunctionToLocalFunctionAdapter<T> >
    {
    public:
      typedef typename T::Traits Traits;
 
      GridFunctionToLocalFunctionAdapter (const T& t_,
                                          const typename Traits::ElementType& e_)
        : t(t_), e(e_) {}
 
      inline void evaluate (const typename Traits::DomainType& x,
                            typename Traits::RangeType& y) const
      {
        t.evaluate(e,x,y);
      }
 
    private:
      const T& t;
      const typename Traits::ElementType& e;
    };
 
 
    template<class T>
    class SelectComponentAdapter : public FunctionInterface<FunctionTraits<typename T::Traits::DomainFieldType,T::Traits::dimDomain,typename T::Traits::DomainType,typename T::Traits::RangeFieldType,1,Dune::FieldVector<typename T::Traits::RangeFieldType,1> > , SelectComponentAdapter<T> >
    {
      typedef FunctionInterface<FunctionTraits<typename T::Traits::DomainFieldType,T::Traits::dimDomain,typename T::Traits::DomainType,typename T::Traits::RangeFieldType,1,Dune::FieldVector<typename T::Traits::RangeFieldType,1> > , SelectComponentAdapter<T> > BaseT;
    public:
      typedef typename BaseT::Traits Traits;
 
      SelectComponentAdapter (const T& t_, int k_) : t(t_), k(k_) {}
 
      inline void evaluate (const typename Traits::DomainType& x,
                            typename Traits::RangeType& y) const
      {
        typename T::Traits::RangeType Y;
        t.evaluate(x,Y);
        y = Y[k];
      }
 
      void select (int k_)
      {
        k = k_;
      }
 
    private:
      const T& t;
      int k;
    };
 
    template<class T>
    class BoundaryGridFunctionSelectComponentAdapter
      : public BoundaryGridFunctionInterface<BoundaryGridFunctionTraits<typename T::Traits::GridViewType,
                                                                        typename T::Traits::RangeFieldType,1,
                                                                        Dune::FieldVector<typename T::Traits::RangeFieldType,1> > ,
                                             BoundaryGridFunctionSelectComponentAdapter<T> >
    {
      typedef BoundaryGridFunctionInterface<BoundaryGridFunctionTraits<typename T::Traits::GridViewType,
                                                                       typename T::Traits::RangeFieldType,1,
                                                                       Dune::FieldVector<typename T::Traits::RangeFieldType,1> > ,
                                            BoundaryGridFunctionSelectComponentAdapter<T> > BaseT;
    public:
      typedef typename BaseT::Traits Traits;
 
      BoundaryGridFunctionSelectComponentAdapter (const T& t_, int k_) : t(t_), k(k_) {}
 
      template<typename I>
      inline void evaluate (const IntersectionGeometry<I>& ig,
                            const typename Traits::DomainType& x,
                            typename Traits::RangeType& y) const
      {
        typename T::Traits::RangeType Y;
        t.evaluate(ig,x,Y);
        y = Y[k];
      }
 
      inline const typename Traits::GridViewType& getGridView () const
      {
        return t.getGridView();
      }
 
 
      void select (int k_)
      {
        k = k_;
      }
 
    private:
      const T& t;
      int k;
    };
 
 
 
  } // namespace PDELab
} // namespace Dune
 
#endif