42 #include <boost/math/special_functions/gamma.hpp>
46 namespace LibUtilities
54 if (lhsPointsKey < rhsPointsKey)
58 if (lhsPointsKey != rhsPointsKey)
104 boost::shared_ptr<Basis> returnval(
new Basis(bkey));
105 returnval->Initialize();
227 D = &(
m_points->GetD()->GetPtr())[0];
243 for (p=0; p<numModes; ++p, mode += numPoints)
247 scal = sqrt(0.5*(2.0*p+1.0));
248 for(i = 0; i < numPoints; ++i)
254 Blas::Dgemm(
'n',
'n',numPoints,numModes,numPoints,1.0,D,numPoints,
272 for(
int p = 0; p < numModes; ++p )
274 for(
int q = 0; q < numModes - p; ++q, mode += numPoints )
277 for(
int j = 0; j < numPoints; ++j )
279 mode[j] *= sqrt(p+q+1.0)*pow(0.5*(1.0 - z[j]), p);
285 Blas::Dgemm(
'n',
'n',numPoints,numModes*(numModes+1)/2,numPoints,1.0,D,numPoints,
302 int P = numModes - 1, Q = numModes - 1, R = numModes - 1;
305 for(
int p = 0; p <= P; ++p )
307 for(
int q = 0; q <= Q - p; ++q )
309 for(
int r = 0; r <= R - p - q; ++r, mode += numPoints )
312 for(
int k = 0; k < numPoints; ++k )
315 mode[k] *= pow(0.5*(1.0 - z[k]), p+q);
318 mode[k] *= sqrt(r+p+q+1.5);
325 Blas::Dgemm(
'n',
'n',numPoints,numModes*(numModes+1)*
326 (numModes+2)/6,numPoints,1.0, D, numPoints,
338 for(i = 0; i < numPoints; ++i)
341 m_bdata[numPoints + i] = 0.5*(1+z[i]);
344 mode =
m_bdata.data() + 2*numPoints;
346 for(p = 2; p < numModes; ++p, mode += numPoints)
350 for(i = 0; i < numPoints; ++i)
357 Blas::Dgemm(
'n',
'n',numPoints,numModes,numPoints,1.0,D,
384 for(i = 0; i < numPoints; ++i)
386 m_bdata[0*numPoints + i] = 0.5*(1-z[i]);
387 m_bdata[1*numPoints + i] = 0.5*(1+z[i]);
390 mode =
m_bdata.data() + 2*numPoints;
392 for(q = 2; q < numModes; ++q, mode+=numPoints)
396 for(i = 0; i < numPoints; ++i)
403 for(i = 0; i < numPoints; ++i)
405 mode[i] = 0.5*(1-z[i]);
410 for(q = 2; q < numModes; ++q, mode+=numPoints)
414 for(i = 0; i < numPoints; ++i)
422 one_p_z =
m_bdata.data()+numPoints;
424 for(p = 2; p < numModes; ++p)
426 for(i = 0; i < numPoints; ++i)
428 mode[i] =
m_bdata[i]*one_m_z_pow[i];
434 for(q = 1; q < numModes-p; ++q, mode+=numPoints)
438 for(i = 0; i < numPoints; ++i)
440 mode[i] *= one_m_z_pow[i]*one_p_z[i];
445 Blas::Dgemm(
'n',
'n',numPoints,numModes*(numModes+1)/2,
446 numPoints,1.0,D,numPoints,
491 for(p = 0; p < numModes; ++p)
493 N = numPoints*(numModes-p)*(numModes-p+1)/2;
495 B_offset += numPoints*(numModes-p);
500 Blas::Dgemm(
'n',
'n',numPoints,
501 numModes*(numModes+1)*(numModes+2)/6,
502 numPoints,1.0,D,numPoints,
514 for (p=0; p<numModes; ++p, mode += numPoints)
516 for(q = 0; q < numPoints; ++q)
518 mode[q] =
Polylib::hglj(p, z[q], zp.data(), numModes, 0.0, 0.0);
523 Blas::Dgemm(
'n',
'n', numPoints, numModes, numPoints, 1.0,
524 D, numPoints,
m_bdata.data(), numPoints, 0.0,
535 for (p=0; p<numModes; ++p,mode += numPoints)
537 for(q = 0; q < numPoints; ++q)
539 mode[q] =
Polylib::hgj(p, z[q], zp.data(), numModes, 0.0, 0.0);
544 Blas::Dgemm(
'n',
'n', numPoints, numModes, numPoints, 1.0,
545 D, numPoints,
m_bdata.data(), numPoints, 0.0,
552 ASSERTL0(numModes%2==0,
"Fourier modes should be a factor of 2");
554 for(i = 0; i < numPoints; ++i)
562 for (p=1; p < numModes/2; ++p)
564 for(i = 0; i < numPoints; ++i)
566 m_bdata[ 2*p *numPoints+i] = cos(p*M_PI* (z[i]+1) );
567 m_bdata[(2*p+1)*numPoints+i] = -sin(p*M_PI* (z[i]+1) );
569 m_dbdata[ 2*p *numPoints+i] = -p*M_PI*sin(p*M_PI* (z[i]+1) );
570 m_dbdata[(2*p+1)*numPoints+i] = -p*M_PI*cos(p*M_PI* (z[i]+1) );
580 for(i = 0; i < numPoints; ++i)
582 m_bdata[i] = cos(M_PI* (z[i]+1) );
583 m_bdata[numPoints+i] = -sin(M_PI* (z[i]+1) );
585 m_dbdata[i] = -M_PI*sin(M_PI* (z[i]+1) );
586 m_dbdata[numPoints+i] = -M_PI*cos(M_PI* (z[i]+1) );
589 for (p=1; p < numModes/2; ++p)
591 for(i = 0; i < numPoints; ++i)
593 m_bdata[ 2*p *numPoints+i] = 0.;
594 m_bdata[(2*p+1)*numPoints+i] = 0.;
597 m_dbdata[(2*p+1)*numPoints+i] = 0.;
618 for (p=0,scal = 1; p<numModes; ++p,mode += numPoints)
622 for(i = 0; i < numPoints; ++i)
627 scal *= 4*(p+1)*(p+1)/(2*p+2)/(2*p+1);
631 Blas::Dgemm(
'n',
'n', numPoints, numModes, numPoints, 1.0,
632 D, numPoints,
m_bdata.data(), numPoints, 0.0,
639 int P = numModes - 1;
642 for(
int p = 0; p <= P; ++p, mode += numPoints )
644 for(
int i = 0; i < numPoints; ++i )
646 mode[i] = pow(z[i], p);
651 Blas::Dgemm(
'n',
'n', numPoints, numModes, numPoints, 1.0,
652 D, numPoints,
m_bdata.data(), numPoints, 0.0,
657 ASSERTL0(
false,
"Basis Type not known or "
658 "not implemented at this time.");
667 bool returnval =
false;
1D Gauss-Radau-Kronrod-Legendre quadrature points, pinned at x=-1
#define ASSERTL0(condition, msg)
BasisKey m_basisKey
Basis specification.
Principle Modified Functions .
const char *const BasisTypeMap[]
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
double hgj(const int i, const double z, const double *zgj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through.
BasisType GetBasisType() const
Return type of expansion basis.
1D Gauss-Radau-Legendre quadrature points, pinned at x=1
Array< OneD, NekDouble > m_bdata
Basis definition.
Principle Modified Functions .
Lagrange Polynomials using the Gauss points .
PointsType GetPointsType() const
Return type of quadrature.
bool operator>(const BasisKey &lhs, const BasisKey &rhs)
BasisType m_basistype
Expansion type.
1D Gauss-Kronrod-Legendre quadrature points
bool operator!=(const BasisKey &x, const BasisKey &y)
1D Gauss-Gauss-Legendre quadrature points
bool operator==(const BasisKey &x, const BasisKey &y)
Principle Orthogonal Functions .
BasisManagerT & BasisManager(void)
std::ostream & operator<<(std::ostream &os, const BasisKey &rhs)
Array< OneD, NekDouble > m_dbdata
Derivative Basis definition.
1D Lobatto Kronrod quadrature points
int GetNumModes() const
Return order of basis from the basis specification.
virtual void Initialize()
Fourier Modified expansions with just the real part of the first mode .
int GetTotNumPoints() const
Return total number of points from the basis specification.
Principle Modified Functions .
void GenBasis()
Generate appropriate basis and their derivatives.
PointsSharedPtr m_points
Set of points.
int GetNumPoints() const
Return points order at which basis is defined.
int m_nummodes
Expansion order.
Principle Orthogonal Functions .
PointsManagerT & PointsManager(void)
int GetTotNumModes() const
Principle Orthogonal Functions .
Defines a specification for a set of points.
boost::shared_ptr< NekMatrix< NekDouble > > CalculateInterpMatrix(const BasisKey &tbasis0)
Calculate the interpolation Matrix for coefficient from one base (m_basisKey) to another (tbasis0) ...
bool operator<(const BasisKey &lhs, const BasisKey &rhs)
Fourier Modified expansions with just the imaginary part of the first mode .
PointsKey GetPointsKey() const
Return distribution of points.
Basis()
Private default constructor.
Fourier ModifiedExpansion with just the first mode .
int GetNumPoints() const
Return the number of points from the basis specification.
Legendre Polynomials . Same as Ortho_A.
double hglj(const int i, const double z, const double *zglj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the.
int GetNumModes() const
Returns the order of the basis.
NekManager< BasisKey, NekMatrix< NekDouble >, BasisKey::opLess > m_InterpManager
boost::shared_ptr< Basis > BasisSharedPtr
bool operator()(const BasisKey &lhs, const BasisKey &rhs) const
static boost::shared_ptr< Basis > Create(const BasisKey &bkey)
Returns a new instance of a Basis with given BasisKey.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
int GetTotNumPoints() const
Describes the specification for a Basis.
bool ExactIprodInt() const
Determine if basis has exact integration for inner product.
BasisType GetBasisType() const
Return the type of expansion basis.
void jacobfd(const int np, const double *z, double *poly_in, double *polyd, const int n, const double alpha, const double beta)
Routine to calculate Jacobi polynomials, , and their first derivative, .
bool Collocation() const
Determine if basis has collocation properties.
1D Gauss-Lobatto-Legendre quadrature points
1D Gauss-Radau-Legendre quadrature points, pinned at x=-1