ROL
NodalBasis.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Rapid Optimization Library (ROL) Package
4//
5// Copyright 2014 NTESS and the ROL contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#include<iostream>
11#include"ROL_StdVector.hpp"
13#include"Lagrange.hpp"
14#include"LinearAlgebra.hpp"
15
16#ifndef __NODAL_BASIS__
17#define __NODAL_BASIS__
18
19
20template<class Real>
22
24 ROL::Ptr<Teuchos::LAPACK<int,Real> > lapack_;
25
27 const int ni_;
28
30 const int nq_;
31
32 NodalBasis(ROL::Ptr<Teuchos::LAPACK<int,Real> > lapack, const int ni, const int nq);
34
36 std::vector<Real> xi_;
37
39 std::vector<Real> xq_;
40
42 std::vector<Real> wq_;
43
45 std::vector<Real> L_;
46 std::vector<Real> Lp_;
47
49 ROL::Ptr<Lagrange<Real> > lagrange_;
50
51
52};
53
56template<class Real>
57NodalBasis<Real>::NodalBasis(ROL::Ptr<Teuchos::LAPACK<int,Real> > const lapack, const int ni, const int nq):
58 lapack_(lapack), ni_(ni), nq_(nq), xi_(ni_,0), xq_(nq_,0), wq_(nq_,0), L_(ni_*nq_,0), Lp_(ni_*nq_,0)
59 {
60
61 // Generate Legendre-Gauss-Lobatto nodes and weights (unused)
62 std::vector<Real> ai(ni_,0);
63 std::vector<Real> bi(ni_,0);
64 std::vector<Real> wi(ni_,0);
65 rec_jacobi(this->lapack_,0,0,ai,bi);
66 rec_lobatto(this->lapack_,-1.0,1.0,ai,bi);
67 gauss(this->lapack_,ai,bi,xi_,wi);
68
69 // Generate Legendre-Gauss nodes and weights
70 std::vector<Real> aq(nq_,0);
71 std::vector<Real> bq(nq_,0);
72 rec_jacobi(this->lapack_,0,0,aq,bq);
73 gauss(this->lapack_,aq,bq,xq_,wq_);
74
75 std::vector<Real> e(ni_,0);
76 std::vector<Real> ell(nq,0);
77
78 lagrange_ = ROL::makePtr<Lagrange<Real>>(xi_,xq_);
79
80 // Loop over canonical vectors
81 for(int i=0;i<ni_;++i) {
82
83 lagrange_->interpolant(i,ell);
84
85 std::copy(ell.begin(),ell.end(),L_.begin()+i*nq_);
86
87
88 lagrange_->derivative(i,ell);
89 std::copy(ell.begin(),ell.end(),Lp_.begin()+i*nq_);
90
91 // Rezero the vector
92 std::fill(e.begin(),e.end(),0);
93 }
94
95}
96
97template<class Real>
100
101#endif
void rec_jacobi(ROL::Ptr< Teuchos::LAPACK< int, Real > > lapack, const double alpha, const double beta, std::vector< Real > &a, std::vector< Real > &b)
Generate the Jacobi polynomial recursion coeffcients .
void gauss(ROL::Ptr< Teuchos::LAPACK< int, Real > > lapack, const std::vector< Real > &a, const std::vector< Real > &b, std::vector< Real > &x, std::vector< Real > &w)
Compute the Gauss quadrature nodes and weights for the polynomials generated by the recurrence coeffi...
void rec_lobatto(ROL::Ptr< Teuchos::LAPACK< int, Real > > const lapack, const double xl1, const double xl2, std::vector< Real > &a, std::vector< Real > &b)
Modify the given recurrence coefficients so that the set of zeros of the maximal order polynomial inc...
std::vector< Real > wq_
std::vector< Real > Lp_
NodalBasis(ROL::Ptr< Teuchos::LAPACK< int, Real > > lapack, const int ni, const int nq)
Set up quantities we will need repeatedly.
const int nq_
std::vector< Real > xi_
ROL::Ptr< Lagrange< Real > > lagrange_
Object for working with Lagrange polynomials and their derivatives.
ROL::Ptr< Teuchos::LAPACK< int, Real > > lapack_
const int ni_
std::vector< Real > xq_
std::vector< Real > L_