NodalBasis.hpp

Go to the documentation of this file.

#include<iostream>
#include"ROL_StdVector.hpp"
#include"OrthogonalPolynomials.hpp"
#include"Lagrange.hpp"
#include"LinearAlgebra.hpp"
 
#ifndef __NODAL_BASIS__
#define __NODAL_BASIS__
 
 
template<class Real>
struct NodalBasis{
 
        ROL::Ptr<Teuchos::LAPACK<int,Real> > lapack_;
 
        const int ni_;
 
        const int nq_;
 
        NodalBasis(ROL::Ptr<Teuchos::LAPACK<int,Real> > lapack, const int ni, const int nq);    
        ~NodalBasis();
 
        std::vector<Real> xi_;
 
        std::vector<Real> xq_;
 
        std::vector<Real> wq_;
 
        std::vector<Real> L_;
        std::vector<Real> Lp_;
        
        ROL::Ptr<Lagrange<Real> > lagrange_; 
 
         
};
 
template<class Real>
NodalBasis<Real>::NodalBasis(ROL::Ptr<Teuchos::LAPACK<int,Real> > const lapack, const int ni, const int nq): 
    lapack_(lapack), ni_(ni), nq_(nq), xi_(ni_,0), xq_(nq_,0), wq_(nq_,0), L_(ni_*nq_,0), Lp_(ni_*nq_,0)
    {
 
    // Generate Legendre-Gauss-Lobatto nodes and weights (unused)
    std::vector<Real> ai(ni_,0);
    std::vector<Real> bi(ni_,0);
    std::vector<Real> wi(ni_,0);
    rec_jacobi(this->lapack_,0,0,ai,bi); 
    rec_lobatto(this->lapack_,-1.0,1.0,ai,bi);
    gauss(this->lapack_,ai,bi,xi_,wi);
 
    // Generate Legendre-Gauss nodes and weights
    std::vector<Real> aq(nq_,0);
    std::vector<Real> bq(nq_,0);
    rec_jacobi(this->lapack_,0,0,aq,bq); 
    gauss(this->lapack_,aq,bq,xq_,wq_);
 
    std::vector<Real> e(ni_,0);
    std::vector<Real> ell(nq,0);
 
    lagrange_ = ROL::makePtr<Lagrange<Real>>(xi_,xq_);
 
    // Loop over canonical vectors
    for(int i=0;i<ni_;++i) {
 
        lagrange_->interpolant(i,ell);
 
        std::copy(ell.begin(),ell.end(),L_.begin()+i*nq_);
      
 
        lagrange_->derivative(i,ell);
        std::copy(ell.begin(),ell.end(),Lp_.begin()+i*nq_);
 
        // Rezero the vector 
        std::fill(e.begin(),e.end(),0);
    }     
 
}
 
template<class Real>
NodalBasis<Real>::~NodalBasis(){
}
 
#endif