gausspoints.hh

#ifndef DUNE_FEM_GAUSSPOINTS_HH
#define DUNE_FEM_GAUSSPOINTS_HH
 
#include <cassert>
#include <vector>
 
#include <dune/common/visibility.hh>
 
namespace Dune
{
 
  namespace Fem
  {
 
    class QuadPtsBase
    {
    protected:
      std::vector< std::vector< double > > G; //[MAXP+1][MAXP]; // positions of Gauss points
      std::vector< std::vector< double > > W; //[MAXP+1][MAXP]; // weights associated with points
      std::vector< int > O;                   //[MAXP+1];       // order of the rule
 
      QuadPtsBase ( const int maxp )
       : G( maxp+1, std::vector<double>(maxp,0.0)),
         W( maxp+1, std::vector<double>(maxp,0.0)),
         O( maxp+1, 0 )
      {}
 
    public:
      double point ( int m, int i ) const
      {
        assert(m > 0 && i < m);
        return G[m][i];
      }
 
      double weight ( int m, int i ) const
      {
        assert(m > 0 && i < m);
        return W[m][i];
      }
 
      int order ( int m ) const
      {
        return O[m];
      }
 
      int power ( int y, int d ) const
      {
        int m = 1;
        for( int i = 0; i < d; ++i )
          m *= y;
        return m;
      }
    };
 
 
    class GaussPts : public QuadPtsBase
    {
    public:
      static const int MAXP=10;
 
      static const int highestOrder=19;
 
    protected:
      using QuadPtsBase :: G;
      using QuadPtsBase :: W;
      using QuadPtsBase :: O;
 
    public:
      GaussPts ();
    };
 
    class ModifiedNewtonCotes : public QuadPtsBase
    {
    public:
      static const int MAXP=10;
 
      static const int highestOrder=10;
 
    protected:
      using QuadPtsBase :: G;
      using QuadPtsBase :: W;
      using QuadPtsBase :: O;
 
    public:
      ModifiedNewtonCotes ();
    };
  } // namespace Fem
 
} // namespace Dune
 
#include "gausspoints_implementation.hh"
#include "modifiednewtoncotes_implementation.hh"
#endif // #ifndef DUNE_FEM_GAUSSPOINTS_HH