tpzgaussrule.h

Go to the documentation of this file.

#ifndef TPZGAUSSRULE_H
#define TPZGAUSSRULE_H

#include "pzvec.h"

#define PZINTEGRAL_MAXITERATIONS_ALLOWED 100000

class TPZGaussRule {
                int fType;

                friend class TPZIntRuleP3D;

    friend class TPZIntRuleList;
                
    int            fNumInt;
    TPZVec<long double>         fLocation;
    TPZVec<long double>         fWeight;

                long double fAlpha;
                long double fBeta;
    
    int fOrder;

    TPZGaussRule(int order,int type = 0,long double alpha = 0.0L,long double beta = 0.0L);
    
    TPZGaussRule(const TPZGaussRule &copy);
    
    TPZGaussRule &operator=(const TPZGaussRule &copy);
    
    ~TPZGaussRule();

public:
    enum {NRULESLOBATTO_ORDER = 645, NRULESLEGENDRE_ORDER = 1000};
                
                int NInt() const{ return fNumInt;}
                
                int Type() { return fType; }
                
                long double Loc(int i) const;
                long double W(int i) const;

                long double Alpha() { return fAlpha; }
                long double Beta() { return fBeta; }
                
    int Order()
    {
        return fOrder;
    }
                void SetParametersJacobi(long double alpha, long double beta) {
                                fAlpha = alpha; fBeta = beta;
                }

                void SetType(int &type,int order);

                void Print(std::ostream & out = std::cout);
                
protected:
                int ComputingGaussLegendreQuadrature(int order);
                void ComputingGaussLegendreQuadrature(int *npoints,TPZVec<long double> &Location,TPZVec<long double> &Weight);
                int ComputingGaussLobattoQuadrature(int order);
//              /** 
//               * @brief Computes the points and weights for Gauss Lobatto Quadrature over the parametric 1D element \f$ [-1.0, 1.0] \f$
//               * It is to integrate functions \f$ f(x) \f$, but one of integration point is \f$-1.0\f$ or \f$1.0\f$ 
//               * @param order Order of the polinomial will be integrated exactly with this cubature rule
//               */
//              int ComputingGaussRadauQuadrature(int order) {
//        DebugStop();
//              }

                int ComputingGaussJacobiQuadrature(int order,long double alpha,long double beta);
                void ComputingGaussJacobiQuadrature(int *npoints,long double alpha, long double beta,TPZVec<long double> &Location,TPZVec<long double> &Weight);

                long double JacobiPolinomial(long double x,int order,long double alpha,long double beta,
                                                                                                                                 TPZVec<long double> &b,TPZVec<long double> &c,long double *dp2,long double *p1);
                
                long double JacobiPolinomial(long double x,int alpha,int beta,unsigned int n);

                int ComputingGaussChebyshevQuadrature(int order);
                
};

long double machinePrecision();

long double gamma(unsigned int n);

long double gamma(long double x);

#endif

//              bool CheckCubatureRule();

/*bool TPZGaussRule::CheckCubatureRule() {
 int i;
 TPZVec<REAL> point(3,0.0L);
 long double sum = 0.0L;
 
 // The cubature rule has not zero integration points
 if(!fNumInt)
 return false;
 
 // Checking all the integration points belong at the master element
 for(i=0;i<fNumInt;i++) {
 point[0] = fLocation[i];
 if(!pztopology::TPZLine::IsInParametricDomain(point))
 break;
 sum += fWeight[i];
 }
 // Checking sum of the weights is equal to measure of the master element
 if(i==fNumInt) {
 if(IsZero((REAL)(sum) - pztopology::TPZLine::RefElVolume()))
 return true;
 }
 return false;   // because any integration point is outside of the master element
 }*/