• 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

dgfunction.hh

Go to the documentation of this file.

#ifndef DG_FUNCTION_HH
#define DG_FUNCTION_HH

#include <dune/istl/bvector.hh>
#include <dune/grid/common/quadraturerules.hh>
#include <dune/disc/functions/functions.hh>
#include <dune/disc/shapefunctions/dgspace/monomialshapefunctions.hh>
#include <dune/disc/shapefunctions/dgspace/orthonormalshapefunctions.hh>

namespace Dune {



  template<class GV, class RT, int o,
           template<class,class,int,int> class SFSC = MonomialShapeFunctionSetContainer>
  class DGFunction :
    virtual public GridFunction<typename GV::Grid,RT,1>,
    virtual public GridFunctionDefault<typename GV::Grid,RT,1>,
    virtual public GridFunctionGlobalEvalDefault<GV,RT,1>,
    virtual public FunctionDefault<typename GV::Grid::ctype,RT,GV::Grid::dimension,1>,
    virtual public H1Function<typename GV::Grid::ctype,RT,GV::Grid::dimension,1>,
    virtual public GridFunctionGlobalDifferentiationDefault<GV,RT,1>,
    virtual public L2Function<typename GV::Grid::ctype,RT,GV::Grid::dimension,1>
  {
    typedef typename GV::Grid G;
    typedef typename GV::IndexSet IS;

  protected:
        typedef typename G::ctype DT;

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

        enum {n=G::dimension};

        DGFunction (const DGFunction&);

  public:
    typedef SFSC<DT,RT,G::dimension,o> ShapeFunctionSetContainer;
    typedef typename ShapeFunctionSetContainer::ShapeFunctionSet ShapeFunctionSet;
  protected:
    static const int BlockSize = ShapeFunctionSetContainer::maxsize;
  public:
    typedef RT ResultType;
        typedef DT DomainType;
    typedef
    FieldVector<ResultType,BlockSize> BlockType;
        typedef BlockVector< BlockType > RepresentationType;

    DGFunction (const GV& gv)
      : GridFunctionGlobalEvalDefault<GV,RT,1>(gv)
      , GridFunctionGlobalDifferentiationDefault<GV,RT,1>(gv)
      , gv_(gv)
    {
      coeff_.resize(gv.indexSet().size(0), false);
    }

        ~DGFunction () {}


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

        RT value = 0;
        const ShapeFunctionSet & s = shapefnkts_(e.type(),o);
        int eid = gv_.indexSet().index(e);
        for (int i=0; i<BlockSize; ++i)
        {
          value += coeff_[eid][i] * s[i].evaluateFunction(0, xi);
        }
        return value;
      }

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

        const ShapeFunctionSet & s = shapefnkts_(e.type(),o);

        typedef Dune::FieldVector<DT,n> Gradient;
        typedef Dune::FieldMatrix<DT,n,n> JacobianInverseMatrix;
        Gradient grad(0);
        Gradient grad_phi_ei;
        JacobianInverseMatrix invj;
        invj = e.geometry().jacobianInverseTransposed(xi);
        int eid = gv_.indexSet().index(e);
      
        for (int i=0; i<BlockSize; ++i)
        {
          Gradient tmp(0.0);
          for (int d=0; d<n; d++)
            tmp[d] = s[i].evaluateDerivative(0, d, xi);
          grad_phi_ei = 0;
          invj.umv(tmp,grad_phi_ei);
          for (int d=0; d<n; d++)
            grad[d] += grad_phi_ei[d] * coeff_[eid][i];
        }

        for (int u=0; u<n; u++)
          if (d[u] == 1) return grad[u];
        DUNE_THROW(NotImplemented,"Can not evaluate this derivative: d=" << d);
      }


        void interpolate (const C0Function<DT,RT,n,1>& u)
      {
        typedef typename GV::template Codim<0>::Iterator Iterator;
        
        Iterator it = gv_.template begin<0> ();
        const Iterator endit = gv_.template end<0> ();
    
        for(; it != endit; ++it)
        {
          FieldVector<RT,BlockSize> l2_b(0.0);
          FieldMatrix<RT,BlockSize,BlockSize> l2_A(0.0);
        
          const ShapeFunctionSet & s = shapefnkts_(it->type(),o);
          int eid = gv_.IndexSet().index(*it);
          coeff_[eid] = 0;
        
          // Get quadrature rule
          typedef QuadratureRule<DT,n> Quad;
          typedef typename QuadratureRule<DT,n>::const_iterator QuadIterator;
          const Quad& quad =
            QuadratureRules<DT,n>::rule(it->type(),
                                        o*getGeometryOrder(it->type()));
        
          for (QuadIterator q=quad.begin(); q!=quad.end(); ++q)
          {
            FieldVector<DT,n> c = q->position();
            FieldVector<DT,n> cg = it->geometry().global(c);
            for(int i=0; i<BlockSize; i++)
            {
              // b_i = \int\limits_{\Omega} exact * Phi_i dV
              RT det = it->geometry().integrationElement(c);
              RT value = u.eval(0, cg);
              RT phi_i = s[i].evaluateFunction(0, c);
              l2_b[i] += q->weight() * phi_i * value * det;
              for(int j=0; j<BlockSize; j++)
              {
                // A_{i,j} = \int\limits_{\Omega} Phi_i * Phi_j dV
                RT phi_j = s[j].evaluateFunction(0, c);
                l2_A[i][j] += q->weight() * phi_i * phi_j * det;
              }
            }
          }
          // solve localy
          fm_solve(l2_A, coeff_[eid], l2_b);
        }
      }


    const RepresentationType& operator* () const
      {
        return coeff_;
      }


    RepresentationType& operator* ()
      {
        return coeff_;
      }

    void preAdapt ()
      {
        DUNE_THROW(NotImplemented, "preAdapt()");
      }

    void postAdapt ()
      {
        DUNE_THROW(NotImplemented, "postAdapt()");
      }
 
    typename VTKWriter<GV>::VTKFunction *
    vtkFunction (std::string s) const
      {
        return new VTKGridFunctionWrapper<GV,RT,1>(*this,s);
      }
        
  protected:
    int getGeometryOrder(const Dune::GeometryType & t)
      {
        if (t.isLine()) return 1;
        if (t.isSimplex()) return 1;
        if (t.isCube()) return t.dim();
        if (t.isPyramid()) return 2;
        if (t.isPrism()) return 2;
        DUNE_THROW(Dune::Exception,"Unknown GeometryType " << t);
      }
    
    // The GridView.  Don't use a reference here since Grid::levelView returns
    // a temporary
    const GV gv_;

    // reference to local shapefunctionset
    ShapeFunctionSetContainer shapefnkts_;

    // and a dynamically allocated vector
    RepresentationType coeff_;
  };

  template<class G, class RT, int o,
           template<class,class,int,int> class SFSC = MonomialShapeFunctionSetContainer>
  class LeafDGFunction : public DGFunction<typename G::LeafGridView,RT,o,SFSC>
  {
  public:
    LeafDGFunction (const G& grid) :
      GridFunctionGlobalEvalDefault<typename G::LeafGridView,RT,1>(grid.leafView()),
      GridFunctionGlobalDifferentiationDefault<typename G::LeafGridView,RT,1>(grid.leafView()),
      DGFunction<typename G::LeafGridView,RT,o,SFSC>(grid.leafView())
      {}
  };

  template<class G, class RT, int o,
           template<class,class,int,int> class SFSC = MonomialShapeFunctionSetContainer>
  class LevelDGFunction : public DGFunction<typename G::LevelGridView,RT,o,SFSC>
  {
  public:
    LevelDGFunction (const G& grid, int level) :
        GridFunctionGlobalEvalDefault<typename G::LevelGridView,RT,1>(grid.levelView(level)),
        GridFunctionGlobalDifferentiationDefault<typename G::LevelGridView,RT,1>(grid.levelView(level)),
        DGFunction<typename G::LevelGridView,RT,o,SFSC>(grid.levelView(level))
    {}
  };

}

#endif // DG_FUNCTION_HH

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