• 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

monomialshapefunctions.hh

Go to the documentation of this file.

#ifndef DUNE_MONOMIALSHAPEFUNCTIONS_HH
#define DUNE_MONOMIALSHAPEFUNCTIONS_HH

#include <iostream>
#include <vector>

#include <dune/common/fvector.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/stdstreams.hh>
#include <dune/common/misc.hh>
#include <dune/grid/common/grid.hh>

#include <dune/disc/shapefunctions/shapefunctions.hh>

namespace Dune
{

  enum {
    MonomialShapeFunctionDefaultMaxOrder = 5
  };

  template<typename C,typename T, int dim>
  class MonomialShapeFunction;
  template<typename C,typename T, int dim>
  class MonomialShapeFunctionSet;
  template<typename C,typename T, int dim,
           int maxOrder=MonomialShapeFunctionDefaultMaxOrder>
  class MonomialShapeFunctionSetContainer;

  template<typename C,typename T, int d>
  class MonomialShapeFunction : public ShapeFunction<C,T,d,1>
  {
  public:
    enum { dim=d  };
    typedef C CoordType;
    typedef T ResultType;
    
    MonomialShapeFunction(int n, const FieldVector<int, dim> & exp):
      n_(n), exp_(exp) {}
    
    virtual ResultType
    evaluateFunction (int, const FieldVector<CoordType,dim>& x) const
      {
        ResultType phi = 1.0;
        // phi= x^(alpha-beta) * y^beta
        // The compiler should be able to eliminate this loop
        for (int c=0; c<dim; c++)
        {
          phi *= power(x[c],exp_[c]);
        }
        return phi;
      } 

    virtual ResultType
    evaluateDerivative(int, int dir, const FieldVector<CoordType,dim>& x) const
      { 
        ResultType dphi = 1.0;
        // The compiler should be able to eliminate this loop
        for (int c=0; c<dim; c++)
        {
          if (c==dir)
            dphi *= exp_[c] * power(x[c],exp_[c]-1);
          else
            dphi *= power(x[c],exp_[c]);
        }
        return dphi;
      }

    virtual int localindex (int comp) const { return n_; };

    virtual int codim () const { return 0; };

    virtual int entity () const { return 0; };

    virtual int entityindex () const { return 0; };

    void print (std::ostream& s) const
    {
      static const char names[] =
        { 'x', 'y', 'z', 'u', 'v', 'w', 'q', 'r', 's' };
      assert( dim < 9 ); /* avoid buffer overrun un names */;
      bool all_zero = true;
      for (int c=0; c<dim; c++)
      {
        if (exp_[c] > 0) {
          if (exp_[c] == 1) {
            s << names[c] << " ";
          }
          else {
            s << names[c] << "^" << exp_[c] << " ";
          }
          all_zero = false;
        }
      }
      if (all_zero) {
        s << 1 <<" ";
      }
    };
  private:
    int n_;
    FieldVector<int, dim> exp_;
    ResultType power(CoordType x, int p) const {
      if (p <= 0)
        return 1.0;
      return x*power(x, p-1);
    }
  };

  template<typename C,typename T, int d>
  inline std::ostream& operator<<(std::ostream& os,
                                  const MonomialShapeFunction<C,T,d>& s)
  {
    s.print(os);
    return os;
  }
  
  template <typename C,typename T, int dim, int c>
  struct MonomialShapeFunctionSetCreator
  {
    typedef MonomialShapeFunction<C,T,dim> ShapeFunction;
    typedef FieldVector<int, dim> expV;
    typedef std::vector< ShapeFunction* > sfV;
    static void create (int bound, int& cnt, expV& exp, sfV& SFnts)
      {
        for (int a=0; a<=bound; a++)
        {
          exp[c-1] = bound-a;
          MonomialShapeFunctionSetCreator<C,T,dim,c+1>::create(a,cnt,exp,SFnts);
        }
      }
  };
  
  template <typename C,typename T, int dim>
  struct MonomialShapeFunctionSetCreator<C,T,dim,0>
  {
    typedef MonomialShapeFunction<C,T,dim> ShapeFunction;
    typedef FieldVector<int, dim> expV;
    typedef std::vector< ShapeFunction* > sfV;
    static int create (int order, sfV& SFnts)
      {
        expV exp(0);
        int cnt = 0;
        for (int a=0; a<=order; a++)
        {
          MonomialShapeFunctionSetCreator<C,T,dim,1>::create(a,cnt,exp,SFnts);
        }
        return cnt;
      }
  };
  
  template <typename C,typename T, int dim>
  struct MonomialShapeFunctionSetCreator<C,T,dim,dim>
  {
    typedef MonomialShapeFunction<C,T,dim> ShapeFunction;
    typedef FieldVector<int, dim> expV;
    typedef std::vector< ShapeFunction* > sfV;
    static void create (int bound, int& cnt, expV& exp, sfV& SFnts)
      {
          exp[dim-1] = bound;
          SFnts[cnt] = new ShapeFunction(cnt, exp);
          for (int i=0; i<dim; i++)
            dvverb << exp[i] << " ";
          dvverb << std::endl;
          cnt++;
      }
  };
  
  template<typename C,typename T, int d>
  class MonomialShapeFunctionSet : public ShapeFunctionSet<C, T, d, 1>
  {
  public:
    enum { dim=d  };
    typedef C CoordType;
    typedef T ResultType;
    typedef MonomialShapeFunction<C,T,dim> ShapeFunction;

    MonomialShapeFunctionSet(int order):
      order_(order), n_(sz(order)), shapeFunctions(n_)
      {
        dvverb << "Constructing monomial shape function set of order "
               << order << std::endl;
#ifndef  NDEBUG
        size_t i =
#endif
          MonomialShapeFunctionSetCreator<C,T,dim,0>::
          create(order,shapeFunctions);
        assert (n_==i);
      }

    virtual ~MonomialShapeFunctionSet()
      {
        for (size_t i=0; i < n_; i++)
        {
          if (shapeFunctions[i]) delete shapeFunctions[i];
        }
      }

    virtual int size() const
      {
        assert(shapeFunctions.size() == n_);
        return n_;
      }

    virtual int size (int entity, int codim) const
      {
        if (codim == 0) return n_;
        return 0;
      };

    virtual const ShapeFunction& operator[](int i) const
      {
        assert (i < (int)n_);
        assert (n_ == shapeFunctions.size());
        return *shapeFunctions[i];
      }

    virtual int order() const
      {
        return order_;
      }

    virtual GeometryType type () const
      {
        static GeometryType t(GeometryType::cube, d);
        return t;
      };
    
  private:
    const int order_;
    const size_t n_;
    static int sz(int order)
      {
        int s;
        s = 1;
        for (int i=1; i<=dim; i++) {
          s = s * (order + i) / i;
        }
        return s;
      }
    std::vector< ShapeFunction* > shapeFunctions;
  };

  template <int dim, int order>
  struct MonomialShapeFunctionSetSize
  {
    enum
    {
      maxSize =
      MonomialShapeFunctionSetSize<dim-1,order>::maxSize
      * (order+dim) / dim
    };
  };

  template <int order>
  struct MonomialShapeFunctionSetSize<1,order>
  {
    enum { maxSize = order + 1 };
  };

  template<typename C,typename T, int d, int maxOrder>
  class MonomialShapeFunctionSetContainer :
    public ShapeFunctionSetContainer<C,T,d,1,MonomialShapeFunctionSetSize<d,maxOrder>::maxSize >
  {
  public:
    enum { dim = d  };
    enum { maxsize = MonomialShapeFunctionSetSize<d,maxOrder>::maxSize };
    typedef C CoordType; 
    typedef T ResultType;
    typedef MonomialShapeFunctionSet<C,T,dim> ShapeFunctionSet;

    MonomialShapeFunctionSetContainer() :
      shapeFunctionSets(maxsize)
      {
        for (int order=0; order<=maxOrder; order++)
        {
          shapeFunctionSets[order] =
            new ShapeFunctionSet(order);
        }
                for (int order=maxOrder+1; order<maxsize; order++)
          shapeFunctionSets[order] = 0;
      }

    virtual ~MonomialShapeFunctionSetContainer()
      {
        for (int i=0; i < maxsize; i++)
        {
          if (shapeFunctionSets[i]) delete shapeFunctionSets[i];
        }
      }

        virtual const ShapeFunctionSet& operator() (GeometryType type, int order) const
      {
        assert(order <= maxOrder);
        return *shapeFunctionSets[order];
      }
  private:
    std::vector< ShapeFunctionSet* > shapeFunctionSets;
  };

}//end of namespace

  
#endif // DUNE_MONOMIALSHAPEFUNCTIONS_HH

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