• Home
• About DUNE
  • What is DUNE
  • News
  • Gallery
  • Publications
• Download
  • Core Modules (Stable)
    • Release Notes (v1.1)
    • Release Notes (v1.0)
  • Recent Changes
  • External Modules
  • Dune Logo
• Documentation
  • Installation
  • Using DUNE
  • Class Documentation
• Community
  • People
  • Mailinglists
  • FAQ
• Development
  • Developing DUNE
  • Unstable Modules
    • Download via SVN
    • Class Documentation
  • Bug Tracking System
  • Automated Test Environment

functions.hh

Go to the documentation of this file.

// $Id: functions.hh 520 2008-10-07 22:06:21Z sander $

#ifndef DUNE_FUNCTIONS_HH
#define DUNE_FUNCTIONS_HH

#include<iostream>
#include<dune/common/fvector.hh>
#include<dune/common/exceptions.hh>
#include<dune/grid/common/grid.hh>
#include<dune/grid/utility/hierarchicsearch.hh>
#include<dune/grid/io/file/vtk/vtkwriter.hh>


namespace Dune
{

  const int InfinitelyDifferentiable=32767;


  template<class DT, class RT, int n, int m>
  class FunctionBase
  {
  public:
        typedef DT DomainFieldType;

        typedef RT RangeFieldType;

        enum { DimDomain = n, DimRange = m};
        

        virtual RT eval (int comp, const Dune::FieldVector<DT,n>& x) const = 0;


        virtual void evalall (const Dune::FieldVector<DT,n>& x, Dune::FieldVector<RT,m>& y) const = 0;

        virtual ~FunctionBase () {}
  };


  template<class DT, class RT, int n, int m>
  class FunctionDefault : virtual public FunctionBase<DT,RT,n,m>
  {
  public:

        virtual void evalall (const Dune::FieldVector<DT,n>& x, Dune::FieldVector<RT,m>& y) const 
        {
          for (int i=0; i<m; ++i)
                y[i] = eval(i,x);
        }
  };






  template<class DT, class RT, int n, int m>
  class DifferentiableFunction : virtual public FunctionBase<DT,RT,n,m>
  {
  public:

        virtual RT derivative (int comp, const Dune::FieldVector<int,n>& d, const Dune::FieldVector<DT,n>& x) const = 0;


        virtual int order () const = 0;
  };



  template<class DT, class RT, int n, int m>
  class DifferentiableFunctionDefault : virtual public DifferentiableFunction<DT,RT,n,m>
  {
  public:

        RT derivative (int comp, const Dune::FieldVector<int,n>& d, const Dune::FieldVector<DT,n>& x) const
        {
          for (int i=0; i<n; ++i)
                if (d[i]>0)
                  {
                        DT eps = 1E-5; // should be type dependent
                        if (sizeof(DT)>4) eps = 1E-10;
                        DT delta = eps*std::abs(x[i])+eps;
                        Dune::FieldVector<int,n> dd(d); dd[i] = dd[i]-1;
                        Dune::FieldVector<DT,n> xx(x); xx[i] += delta;
                        RT up = derivative(comp,dd,xx);
                        xx[i] = x[i]-delta;
                        RT down = derivative(comp,dd,xx);
                        return (up-down)/(2*delta);
                  }
          // d contains all zeros, evaluate function
          return eval(comp,x);
        }
  };




  template<class G, class RT, int m>
  class GridFunction : virtual public FunctionBase<typename G::ctype,RT,G::dimension,m>
  {
        typedef typename G::ctype DT;

        enum {n=G::dimension};

        typedef typename G::Traits::template Codim<0>::Entity Entity;

  public:

        virtual RT evallocal (int comp, const Entity& e, const Dune::FieldVector<DT,n>& xi) const = 0;


        virtual void evalalllocal (const Entity& e, const Dune::FieldVector<DT,G::dimension>& xi, 
                                                           Dune::FieldVector<RT,m>& y) const = 0;
  };



  template<class G, class RT, int m>
  class GridFunctionDefault : virtual public GridFunction<G,RT,m>
  {
        typedef typename G::ctype DT;

        enum {n=G::dimension};

        typedef typename G::Traits::template Codim<0>::Entity Entity;

  public:

        RT evallocal (int comp, const Entity& e, const Dune::FieldVector<DT,n>& xi) const
        {
          return eval(comp,e.geometry().global(xi));
        }



        void evalalllocal (const Entity& e, const Dune::FieldVector<DT,G::dimension>& xi, 
                                           Dune::FieldVector<RT,m>& y) const
        {
          for (int i=0; i<m; ++i)
          y[i] = evallocal(i,e,xi);
        }

  };


  template<class GV, class RT, int m>
  class GridFunctionGlobalEvalDefault : virtual public GridFunction<typename GV::Grid,RT,m>
  {
        typedef typename GV::Grid G;

        typedef typename G::ctype DT;

        enum {n=G::dimension};

        typedef typename G::Traits::template Codim<0>::Entity Entity;

        typedef typename G::Traits::template Codim<0>::EntityPointer EntityPointer;

      typedef typename GV::IndexSet IS;
  public:

      GridFunctionGlobalEvalDefault(const GV& _gv) : hsearch(_gv.grid(),_gv.indexSet()) {};
    
        virtual RT eval (int comp, const Dune::FieldVector<DT,n>& xi) const
        {
      EntityPointer e = hsearch.findEntity(xi);
      return evallocal(comp,*e,e->geometry().local(xi));
        }

  protected:
    HierarchicSearch<G,IS> hsearch;
  };


  template<class G, class RT, int m>
  class DifferentiableGridFunction : virtual public GridFunction<G,RT,m>, 
                                                                         virtual public DifferentiableFunction<typename G::ctype,RT,G::dimension,m>
  {
        typedef typename G::ctype DT;

        enum {n=G::dimension};

        typedef typename G::Traits::template Codim<0>::Entity Entity;

  public:

        virtual RT derivativelocal (int comp, const Dune::FieldVector<int,n>& d, 
                                                                const Entity& e, const Dune::FieldVector<DT,n>& xi) const = 0;
  };
  

  template<class G, class RT, int m>
  class DifferentiableGridFunctionDefault : virtual public DifferentiableGridFunction<G,RT,m>
  {
        typedef typename G::ctype DT;

        enum {n=G::dimension};

        typedef typename G::Traits::template Codim<0>::Entity Entity;

  public:

        RT derivativelocal (int comp, const Dune::FieldVector<int,n>& d, 
                                                const Entity& e, const Dune::FieldVector<DT,n>& xi) const
        { 
          return derivative(comp,d,e.geometry().global(xi));
        }
  };
  
  template<class GV, class RT, int m>
  class GridFunctionGlobalDifferentiationDefault :
    virtual public DifferentiableGridFunction<typename GV::Grid,RT,m>,
    virtual public GridFunctionGlobalEvalDefault<GV,RT,m>
  {
        typedef typename GV::Grid G;

        typedef typename G::ctype DT;

        enum {n=G::dimension};

        typedef typename G::Traits::template Codim<0>::Entity Entity;

        typedef typename G::Traits::template Codim<0>::EntityPointer EntityPointer;

  public:

    GridFunctionGlobalDifferentiationDefault(const GV & _gv) :
      GridFunctionGlobalEvalDefault<GV,RT,m>(_gv) {};
    
    virtual RT derivative (int comp, const Dune::FieldVector<int, n> &d,
                           const Dune::FieldVector<DT, n> &x) const
        {
      EntityPointer e = this->hsearch.findEntity(x);
      return derivativelocal(comp,d,*e,e->geometry().local(x));
        }
  };


  template<class DT, class RT, int n, int m>
  class C0Function : virtual public FunctionBase<DT,RT,n,m>
  {
  };


  template<class G, class RT, int m>
  class C0GridFunction : virtual public GridFunction<G,RT,m>,
                                                 virtual public C0Function<typename G::ctype,RT,G::dimension,m>
  {
  };
  


  template<class DT, class RT, int n, int m>
  class C1Function : virtual public DifferentiableFunction<DT,RT,n,m>,
                                         virtual public C0Function<DT,RT,n,m>
  {
  public:

        virtual int order () const {return 1;}
  };




  template<class DT, class RT, int n, int m>
  class L2Function : virtual public FunctionBase<DT,RT,n,m>
  {
  };




  template<class DT, class RT, int n, int m>
  class H1Function : virtual public DifferentiableFunction<DT,RT,n,m>,
                                         virtual public L2Function<DT,RT,n,m>
  {
  public:

        virtual int order () const {return 1;}
  };



  template<class DT, class RT, int n, int m>
  class HdivFunction : virtual public DifferentiableFunction<DT,RT,n,m>,
                                         virtual public L2Function<DT,RT,n,m>
  {
  public:

        virtual int order () const {return 1;}
  };




  template<class G, class RT, int m>
  class ElementwiseCInfinityFunction : virtual public DifferentiableGridFunction<G,RT,m>
  {
  public:

        virtual int order () const {return InfinitelyDifferentiable;}
  };
  

  // Wrapper for VTK output of grid functions
  template<class GV, class RT, int m>
  class VTKGridFunctionWrapper : public VTKWriter<GV>::VTKFunction  
  {
      enum {n=GV::dimension};
      typedef typename GV::Grid::ctype DT;
        typedef typename GV::template Codim<0>::Entity Entity;

  public:
        virtual int ncomps () const
        {
          return m;
        }


        virtual double evaluate (int comp, const Entity& e, const Dune::FieldVector<DT,n>& xi) const
        {
          return func.evallocal(comp,e,xi);
        }
          
        // get name
        virtual std::string name () const
        {
          return myname;
        }

        // constructor remembers reference to a grid function
      VTKGridFunctionWrapper (const GridFunction<typename GV::Grid,RT,m>& f, std::string s) : func(f), myname(s)
        {}

        virtual ~VTKGridFunctionWrapper() {}
      
  private:
      const GridFunction<typename GV::Grid,RT,m>& func;
        std::string myname;
  };

  

}
#endif

Generated on 6 Jan 2009 with Doxygen (ver 1.5.1) [logfile].