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/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 comp, const Entity& e,
                             const Dune::FieldVector<ctype,dim>& xi) const
    {
      return v[mapper.map(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
    {
      double min=1E100;
      int imin=-1;
      Dune::GeometryType gt = e.type();
      for (int i=0; i<e.template count<dim>(); ++i)
      {
        Dune::FieldVector<ctype,dim> local =
          Dune::ReferenceElements<ctype,dim>::general(gt)
          .position(i,dim);
        local -= xi;
        if (local.infinity_norm()<min)
        {
          min = local.infinity_norm();
          imin = i;
        }
      }
      return v[mapper.map(e,imin,dim)*ncomps_+mycomp_];
    }
 
    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