[dune-pdelab] Derive from DiscreteGridFunction

[Christoph Grüninger]

Hi,
I want to write the gradient of my solution to a Vtk file. I derived from DiscreteGridFunction and wrote my own evalute function but I can't access the GridFunctionSpace from DiscreteGridFunction because it is private. I could store it in a private field, too. But the derivation would be useless. Where is my mistake?

Cheers
Christoph


My code:

// get gradient from solution
template<typename GFS, typename GCV>
class SolutionVelocity
  : public Dune::PDELab::DiscreteGridFunction<GFS, GCV>
{
  typedef typename Dune::PDELab::GridFunctionInterface<
    Dune::PDELab::GridFunctionTraits<
      typename GFS::Traits::GridViewType,
      typename GFS::Traits::LocalFiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
      GFS::Traits::GridViewType::Grid::dimensionworld,
      Dune::FieldVector<
        typename GFS::Traits::LocalFiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
        GFS::Traits::GridViewType::Grid::dimensionworld> >,
    SolutionVelocity<GFS,GCV> > BaseT;
  typedef typename Dune::PDELab::DiscreteGridFunction<GFS,GCV> DGF;
  
public:
  typedef typename BaseT::Traits Traits;
  typedef typename Traits::GridViewType GV;
  typedef typename Traits::DomainFieldType DF;
  typedef typename Traits::RangeFieldType RF;
  typedef typename Traits::RangeType RangeType;
  typedef typename GFS::Traits::LocalFiniteElementType::
    Traits::LocalBasisType::Traits::JacobianType JacobianType;

  SolutionVelocity (const GFS& gfs, const GCV& gcv)
    : DGF::DiscreteGridFunction(gfs, gcv) {}

  // evaluate
  inline void evaluate (const typename Traits::ElementType& e,
                        const typename Traits::DomainType& x,
                        typename Traits::RangeType& y) const
  {
    const int dim = GV::Grid::dimension;
    typedef K<GV,RF> KType;
    //KType k(DGF::pgfs->gridview());
    KType k(DGF::getGridView());
    DGF::lfs.bind(e); // Hier meckert der Compiler
		DGF::lfs.vread(DGF::xg,DGF::xl);

    // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
    std::vector<JacobianType> js(DGF::lfs.size());
    DGF::lfs.localFiniteElement().localBasis().evaluateJacobian(x,js);
    Dune::FieldVector<DF,dim> gradu(0.0);

    // transform gradient to real element
    const Dune::FieldMatrix<DF,dim,dim> jac = e.geometry().jacobianInverseTransposed(x);
    std::vector<Dune::FieldVector<RF,dim> > gradphi(DGF::lfs.size());
    for (size_t j=0; j<DGF::lfs.size(); j++)
    {
      gradphi[j] = 0.0;
      jac.mv(js[j][0],gradphi[j]);
      gradu.axpy(DGF::xl[j],gradphi[j]);
    }
    typename KType::Traits::RangeType permeability(0.0);
    const Dune::FieldVector<DF,dim> localcenter(0.5);
    k.evaluate(e,localcenter,permeability);
    permeability.mv(gradu,y);
    y *= -1.0;
  }
};