function.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_GRID_IO_FILE_VTK_FUNCTION_HH
#define DUNE_GRID_IO_FILE_VTK_FUNCTION_HH
 
#include <string>
 
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
 
#include <dune/geometry/type.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/multilineargeometry.hh>
 
#include <dune/grid/common/mcmgmapper.hh>
 
namespace Dune
{
 
  //
  //  Base VTKFunction
  //
 
  template< class GridView >
  class VTKFunction
  {
  public:
    typedef typename GridView::ctype ctype;
    enum { dim = GridView::dimension };
    typedef typename GridView::template Codim< 0 >::Entity Entity;
 
    virtual int ncomps () const = 0;
 
 
    virtual double evaluate (int comp, const Entity& e,
                             const Dune::FieldVector<ctype,dim>& xi) const = 0;
 
    virtual std::string name () const = 0;
 
    virtual ~VTKFunction () {}
  };
 
  //
  //  P0VTKFunction
  //
 
 
  template<typename GV, typename V>
  class P0VTKFunction
    : public VTKFunction< GV >
  {
    typedef VTKFunction< GV > Base;
    typedef MultipleCodimMultipleGeomTypeMapper<GV, MCMGElementLayout> Mapper;
 
    const V& v;
    std::string s;
    int ncomps_;
    int mycomp_;
    Mapper mapper;
 
  public:
    typedef typename Base::Entity Entity;
    typedef typename Base::ctype ctype;
    using Base::dim;
 
    virtual int ncomps () const
    {
      return 1;
    }
 
    virtual double evaluate (int, const Entity& e,
                             const Dune::FieldVector<ctype,dim>&) const
    {
      return v[mapper.index(e)*ncomps_+mycomp_];
    }
 
    virtual std::string name () const
    {
      return s;
    }
 
 
    P0VTKFunction(const GV &gv, const V &v_, const std::string &s_,
                  int ncomps=1, int mycomp=0 )
      : v( v_ ),
        s( s_ ),
        ncomps_(ncomps),
        mycomp_(mycomp),
        mapper( gv )
    {
      if (v.size()!=(unsigned int)(mapper.size()*ncomps_))
        DUNE_THROW(IOError, "P0VTKFunction: size mismatch");
    }
 
    virtual ~P0VTKFunction() {}
  };
 
  //
  //  P1VTKFunction
  //
 
 
  template<typename GV, typename V>
  class P1VTKFunction
    : public VTKFunction< GV >
  {
    typedef VTKFunction< GV > Base;
    typedef MultipleCodimMultipleGeomTypeMapper<GV, MCMGVertexLayout> Mapper;
 
    const V& v;
    std::string s;
    int ncomps_;
    int mycomp_;
    Mapper mapper;
 
  public:
    typedef typename Base::Entity Entity;
    typedef typename Base::ctype ctype;
    using Base::dim;
 
    virtual int ncomps () const
    {
      return 1;
    }
 
    virtual double evaluate (int comp, const Entity& e,
                             const Dune::FieldVector<ctype,dim>& xi) const
    {
      const unsigned int dim = Entity::mydimension;
      const unsigned int nVertices = e.subEntities(dim);
 
      std::vector<FieldVector<ctype,1> > cornerValues(nVertices);
      for (unsigned i=0; i<nVertices; ++i)
        cornerValues[i] = v[mapper.subIndex(e,i,dim)*ncomps_+mycomp_];
 
      // (Ab)use the MultiLinearGeometry class to do multi-linear interpolation between scalars
      const MultiLinearGeometry<ctype,dim,1> interpolation(e.type(), cornerValues);
      return interpolation.global(xi);
    }
 
    virtual std::string name () const
    {
      return s;
    }
 
 
    P1VTKFunction(const GV& gv, const V &v_, const std::string &s_,
                  int ncomps=1, int mycomp=0 )
      : v( v_ ),
        s( s_ ),
        ncomps_(ncomps),
        mycomp_(mycomp),
        mapper( gv )
    {
      if (v.size()!=(unsigned int)(mapper.size()*ncomps_))
        DUNE_THROW(IOError,"P1VTKFunction: size mismatch");
    }
 
    virtual ~P1VTKFunction() {}
  };
 
 
} // namespace Dune
 
#endif // DUNE_GRID_IO_FILE_VTK_FUNCTION_HH