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