Represents a basis of a given type. More...

#include <Basis.h>

Public Member Functions
virtual	~Basis ()
	Destructor. More...

int	GetNumModes () const
	Return order of basis from the basis specification. More...

int	GetTotNumModes () const
	Return total number of modes from the basis specification. More...

int	GetNumPoints () const
	Return the number of points from the basis specification. More...

int	GetTotNumPoints () const
	Return total number of points from the basis specification. More...

BasisType	GetBasisType () const
	Return the type of expansion basis. More...

PointsKey	GetPointsKey () const
	Return the points specification for the basis. More...

PointsType	GetPointsType () const
	Return the type of quadrature. More...

const Array< OneD, const NekDouble > &	GetZ () const

const Array< OneD, const NekDouble > &	GetW () const

void	GetZW (Array< OneD, const NekDouble > &z, Array< OneD, const NekDouble > &w) const

const std::shared_ptr< NekMatrix< NekDouble > > &	GetD (Direction dir=xDir) const

const std::shared_ptr< NekMatrix< NekDouble > >	GetI (const Array< OneD, const NekDouble > &x)

const std::shared_ptr< NekMatrix< NekDouble > >	GetI (const BasisKey &bkey)

bool	ExactIprodInt () const
	Determine if basis has exact integration for inner product. More...

bool	Collocation () const
	Determine if basis has collocation properties. More...

const Array< OneD, const NekDouble > &	GetBdata () const
	Return basis definition array m_bdata. More...

const Array< OneD, const NekDouble > &	GetDbdata () const
	Return basis definition array m_dbdata. More...

const BasisKey	GetBasisKey () const
	Returns the specification for the Basis. More...

virtual void	Initialize ()

Static Public Member Functions
static std::shared_ptr< Basis >	Create (const BasisKey &bkey)
	Returns a new instance of a Basis with given BasisKey. More...

Protected Attributes
BasisKey	m_basisKey
	Basis specification. More...

PointsSharedPtr	m_points
	Set of points. More...

Array< OneD, NekDouble >	m_bdata
	Basis definition. More...

Array< OneD, NekDouble >	m_dbdata
	Derivative Basis definition. More...

NekManager< BasisKey, NekMatrix< NekDouble >, BasisKey::opLess >	m_InterpManager

Private Member Functions
	Basis (const BasisKey &bkey)
	Private constructor with BasisKey. More...

	Basis ()
	Private default constructor. More...

std::shared_ptr< NekMatrix< NekDouble > >	CalculateInterpMatrix (const BasisKey &tbasis0)
	Calculate the interpolation Matrix for coefficient from one base (m_basisKey) to another (tbasis0) More...

void	GenBasis ()
	Generate appropriate basis and their derivatives. More...

Static Private Attributes
static bool	initBasisManager []

Detailed Description

Represents a basis of a given type.

Definition at line 211 of file Basis.h.

Constructor & Destructor Documentation

◆ ~Basis()

virtual Nektar::LibUtilities::Basis::~Basis ( )

inlinevirtual

Destructor.

Definition at line 218 of file Basis.h.

219 {

220 };

◆ Basis() [1/2]

Nektar::LibUtilities::Basis::Basis ( const BasisKey & bkey )

private

Private constructor with BasisKey.

Definition at line 95 of file Foundations/Basis.cpp.

References CalculateInterpMatrix(), and m_InterpManager.

                                         :
             m_basisKey(bkey),
             m_points(PointsManager()[bkey.GetPointsKey()]),
             m_bdata(bkey.GetTotNumModes()*bkey.GetTotNumPoints()),
             m_dbdata(bkey.GetTotNumModes()*bkey.GetTotNumPoints())
         {
             m_InterpManager.RegisterGlobalCreator(std::bind(&Basis::CalculateInterpMatrix,this,std::placeholders::_1));
         }

◆ Basis() [2/2]

Nektar::LibUtilities::Basis::Basis ( )

inlineprivate

Private default constructor.

Definition at line 347 of file Basis.h.

References Nektar::ErrorUtil::efatal, Nektar::lhs, LIB_UTILITIES_EXPORT, NEKERROR, Nektar::LibUtilities::BasisKey::operator<, Nektar::LibUtilities::operator<<(), Nektar::LibUtilities::operator>(), and Nektar::rhs.

Referenced by Create().

                    :m_basisKey(NullBasisKey)
             {
                 NEKERROR(ErrorUtil::efatal,
                          "Default Constructor for Basis should not be called");
             }

Member Function Documentation

◆ CalculateInterpMatrix()

std::shared_ptr< NekMatrix< NekDouble > > Nektar::LibUtilities::Basis::CalculateInterpMatrix ( const BasisKey & tbasis0 )

private

Calculate the interpolation Matrix for coefficient from one base (m_basisKey) to another (tbasis0)

Definition at line 124 of file Foundations/Basis.cpp.

References Nektar::LibUtilities::BasisManager(), Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::BasisKey::GetBasisType(), Nektar::LibUtilities::BasisKey::GetNumModes(), and m_basisKey.

Referenced by Basis().

         {
             int dim = m_basisKey.GetNumModes();
             const PointsKey pkey(dim,LibUtilities::eGaussLobattoLegendre);
             BasisKey fbkey(m_basisKey.GetBasisType(),dim,pkey);
             BasisKey tbkey(tbasis0.GetBasisType(),dim,pkey);
 
             // "Constructur" of the basis
             BasisSharedPtr fbasis = BasisManager()[fbkey];
             BasisSharedPtr tbasis = BasisManager()[tbkey];
 
             // Get B Matrices
             Array<OneD, NekDouble> fB_data = fbasis->GetBdata();
             Array<OneD, NekDouble> tB_data = tbasis->GetBdata();
 
             // Convert to a NekMatrix
             NekMatrix<NekDouble> fB(dim,dim,fB_data);
             NekMatrix<NekDouble> tB(dim,dim,tB_data);
 
             // Invert the "to" matrix: tu = tB^(-1)*fB fu = ftB fu
             tB.Invert();
 
             // Compute transformation matrix
             Array<OneD, NekDouble> zero1D(dim*dim,0.0);
             std::shared_ptr< NekMatrix<NekDouble> > ftB(MemoryManager<NekMatrix<NekDouble> >::AllocateSharedPtr(dim,dim,zero1D));
             (*ftB) = tB*fB;
 
             return ftB;
         }

◆ Collocation()

bool Nektar::LibUtilities::Basis::Collocation ( ) const

inline

Determine if basis has collocation properties.

Definition at line 307 of file Basis.h.

             {
                 return m_basisKey.Collocation();
             }

◆ Create()

std::shared_ptr< Basis > Nektar::LibUtilities::Basis::Create ( const BasisKey & bkey )

static

Returns a new instance of a Basis with given BasisKey.

Definition at line 104 of file Foundations/Basis.cpp.

References Basis().

         {
             std::shared_ptr<Basis> returnval(new Basis(bkey));
             returnval->Initialize();
 
             return returnval;
         }

◆ ExactIprodInt()

bool Nektar::LibUtilities::Basis::ExactIprodInt ( void ) const

inline

Determine if basis has exact integration for inner product.

Definition at line 301 of file Basis.h.

             {
                 return m_basisKey.ExactIprodInt();
             }

◆ GenBasis()

void Nektar::LibUtilities::Basis::GenBasis ( )

private

Generate appropriate basis and their derivatives.

The following expansions are generated depending on the enum type defined in m_basisKey.m_basistype:

NOTE: This definition does not follow the order in the Karniadakis & Sherwin book since this leads to a more compact hierarchical pattern for implementation purposes. The order of these modes dictates the ordering of the expansion coefficients.

In the following m_numModes = P

eModified_A:

m_bdata[i + j*m_numpoints] = \( \phi^a_i(z_j) = \left \{ \begin{array}{ll} \left ( \frac{1-z_j}{2}\right ) & i = 0 \\ \\ \left ( \frac{1+z_j}{2}\right ) & i = 1 \\ \\ \left ( \frac{1-z_j}{2}\right )\left ( \frac{1+z_j}{2}\right ) P^{1,1}_{i-2}(z_j) & 2\leq i < P\\ \end{array} \right . \)

eModified_B:

m_bdata[n(i,j) + k*m_numpoints] = \( \phi^b_{ij}(z_k) = \left \{ \begin{array}{lll} \phi^a_j(z_k) & i = 0, & 0\leq j < P \\ \\ \left ( \frac{1-z_k}{2}\right )^{i} & 1 \leq i < P,& j = 0 \\ \\ \left ( \frac{1-z_k}{2}\right )^{i} \left ( \frac{1+z_k}{2}\right ) P^{2i-1,1}_{j-1}(z_k) & 1 \leq i < P,\ & 1\leq j < P-i\ \\ \end{array} \right . , \)

where \( n(i,j) \) is a consecutive ordering of the triangular indices \( 0 \leq i, i+j < P \) where j runs fastest.

eModified_C:

m_bdata[n(i,j,k) + l*m_numpoints] = \( \phi^c_{ij,k}(z_l) = \phi^b_{i+j,k}(z_l) = \left \{ \begin{array}{llll} \phi^b_{j,k}(z_l) & i = 0, & 0\leq j < P & 0\leq k < P-j\\ \\ \left ( \frac{1-z_l}{2}\right )^{i+j} & 1\leq i < P,\ & 0\leq j < P-i,\ & k = 0 \\ \\ \left ( \frac{1-z_l}{2}\right )^{i+j} \left ( \frac{1+z_l}{2}\right ) P^{2i+2j-1,1}_{k-1}(z_k) & 1\leq i < P,& 0\leq j < P-i& 1\leq k < P-i-j \\ \\ \end{array} \right . , \)

where \( n(i,j,k) \) is a consecutive ordering of the triangular indices \( 0 \leq i, i+j, i+j+k < P \) where k runs fastest, then j and finally i.

Orthogonal basis A

\(\tilde \psi_p^a (\eta_1) = L_p(\eta_1) = P_p^{0,0}(\eta_1)\)

Orthogonal basis B

\(\tilde \psi_{pq}^b(\eta_2) = \left ( {1 - \eta_2} \over 2 \right)^p P_q^{2p+1,0}(\eta_2)\) \

Orthogonal basis C

\(\tilde \psi_{pqr}^c = \left ( {1 - \eta_3} \over 2 \right)^{p+q} P_r^{2p+2q+2, 0}(\eta_3)\) \ \

Orthogonal basis C for Pyramid expansion (which is richer than tets)

\(\tilde \psi_{pqr}^c = \left ( {1 - \eta_3} \over 2\right)^{pq} P_r^{2pq+2, 0}(\eta_3)\) \( \mbox{where }pq = max(p+q,0) \)

This orthogonal expansion has modes that are always in the Cartesian space, however the equivalent ModifiedPyr_C has vertex modes that do not lie in this space. If one chooses \(pq = max(p+q-1,0)\) then the expansion will space the same space as the vertices but the order of the expanion in 'r' is reduced by one.

1) Eta_z values are the changing the fastest, then r, q, and finally p. 2) r index increases by the stride of numPoints.

Definition at line 217 of file Foundations/Basis.cpp.

Referenced by Initialize().

         {
             int i,p,q;
             NekDouble scal;
             Array<OneD, NekDouble> modeSharedArray;
             NekDouble *mode;
             Array<OneD, const NekDouble> z;
             Array<OneD, const NekDouble> w;
             const NekDouble *D;
 
             m_points->GetZW(z,w);
 
             D = &(m_points->GetD()->GetPtr())[0];
             int numModes = GetNumModes();
             int numPoints = GetNumPoints();
 
             switch(GetBasisType())
             {
 
             /** \brief Orthogonal basis A
 
             \f$\tilde \psi_p^a (\eta_1) = L_p(\eta_1) = P_p^{0,0}(\eta_1)\f$
 
            */
             case eOrtho_A:
             case eLegendre:
                 mode = m_bdata.data();
             
                 for (p=0; p<numModes; ++p, mode += numPoints)
                 {
                     Polylib::jacobfd(numPoints, z.data(), mode, NULL, p, 0.0, 0.0);
                     // normalise
                     scal = sqrt(0.5*(2.0*p+1.0));
                     for(i = 0; i < numPoints; ++i)
                     {
                         mode[i] *= scal;
                     }
                 }
                 // define derivative basis
                 Blas::Dgemm('n','n',numPoints,numModes,numPoints,1.0,D,numPoints,
                             m_bdata.data(),numPoints,0.0,m_dbdata.data(),numPoints);
                 break;
 
             /** \brief Orthogonal basis B
 
             \f$\tilde \psi_{pq}^b(\eta_2) = \left ( {1 - \eta_2} \over 2 \right)^p P_q^{2p+1,0}(\eta_2)\f$ \\
 
            */
 
             // This is tilde psi_pq in Spen's book, page 105
             // The 3-dimensional array is laid out in memory such that
             // 1) Eta_y values are the changing the fastest, then q and p.
             // 2) q index increases by the stride of numPoints.
             case eOrtho_B:
                 {
                      NekDouble *mode = m_bdata.data();
 
                      for( int p = 0; p < numModes; ++p )
                      {
                          for( int q = 0; q < numModes - p; ++q,  mode += numPoints  )
                          {
                              Polylib::jacobfd(numPoints, z.data(), mode, NULL, q, 2*p + 1.0, 0.0);
                              for( int j = 0; j < numPoints; ++j )
                              {
                                  mode[j] *= sqrt(p+q+1.0)*pow(0.5*(1.0 - z[j]), p);
                              }
                          }
                      }
 
                      // define derivative basis
                      Blas::Dgemm('n','n',numPoints,numModes*(numModes+1)/2,numPoints,1.0,D,numPoints,
                                   m_bdata.data(),numPoints,0.0,m_dbdata.data(),numPoints);
                 }
                 break;
 
             /** \brief Orthogonal basis C
 
                 \f$\tilde \psi_{pqr}^c = \left ( {1 - \eta_3} \over 2 \right)^{p+q} P_r^{2p+2q+2, 0}(\eta_3)\f$ \ \
 
            */
 
             // This is tilde psi_pqr in Spen's book, page 105
             // The 4-dimensional array is laid out in memory such that
             // 1) Eta_z values are the changing the fastest, then r, q, and finally p.
             // 2) r index increases by the stride of numPoints.
             case eOrtho_C:
                 {
                     int P = numModes - 1, Q = numModes - 1, R = numModes - 1;
                     NekDouble *mode = m_bdata.data();
 
                     for( int p = 0; p <= P; ++p )
                     {
                         for( int q = 0; q <= Q - p; ++q )
                         {
                             for( int r = 0; r <= R - p - q; ++r, mode += numPoints )
                             {
                                 Polylib::jacobfd(numPoints, z.data(), mode, NULL, r, 2*p + 2*q + 2.0, 0.0);
                                 for( int k = 0; k < numPoints; ++k )
                                 {
                                     // Note factor of 0.5 is part of normalisation
                                     mode[k] *= pow(0.5*(1.0 - z[k]), p+q);
 
                                     // finish normalisation
                                     mode[k] *= sqrt(r+p+q+1.5);
                                 }
                             }
                         }
                     }
 
                     // Define derivative basis
                     Blas::Dgemm('n','n',numPoints,numModes*(numModes+1)*
                         (numModes+2)/6,numPoints,1.0, D, numPoints,
                         m_bdata.data(),numPoints,0.0,m_dbdata.data(),numPoints);
                 }
                 break;
 
                 /** \brief Orthogonal basis C for Pyramid expansion
                     (which is richer than tets)
                     
                     \f$\tilde \psi_{pqr}^c = \left ( {1 - \eta_3} \over
                     2\right)^{pq} P_r^{2pq+2, 0}(\eta_3)\f$ \f$ \mbox{where
                     }pq = max(p+q,0) \f$
                     
                     This orthogonal expansion has modes that are
                     always in the Cartesian space, however the
                     equivalent ModifiedPyr_C has vertex modes that do
                     not lie in this space. If one chooses \f$pq =
                     max(p+q-1,0)\f$ then the expansion will space the
                     same space as the vertices but the order of the
                     expanion in 'r' is reduced by one.
 
                     1) Eta_z values are the changing the fastest, then
                      r, q, and finally p.  2) r index increases by the
                       stride of numPoints. 
                 */
             case eOrthoPyr_C:
                 {
                     int P = numModes - 1, Q = numModes - 1, R = numModes - 1;
                     NekDouble *mode = m_bdata.data();
 
                     for( int p = 0; p <= P; ++p )
                     {
                         for( int q = 0; q <= Q; ++q )
                         {
                             for( int r = 0; r <= R - std::max(p,q); ++r, mode += numPoints )
                             {
                                 // this offset allows for orthogonal
                                 // expansion to span linear FE space
                                 // of modified basis but means that
                                 // the cartesian polynomial space
                                 // spanned by the expansion is one
                                 // order lower.
                                 //int pq = max(p + q -1,0);
                                 int pq = std::max(p + q,0);
                                 
                                 Polylib::jacobfd(numPoints, z.data(), mode, NULL, r, 2*pq + 2.0, 0.0);
                                 for( int k = 0; k < numPoints; ++k )
                                 {
                                     // Note factor of 0.5 is part of normalisation
                                     mode[k] *= pow(0.5*(1.0 - z[k]), pq);
 
                                     // finish normalisation
                                     mode[k] *= sqrt(r+pq+1.5);
                                 }
                             }
                         }
                     }
 
                     // Define derivative basis
                     Blas::Dgemm('n','n',numPoints,numModes*(numModes+1)*
                         (numModes+2)/6,numPoints,1.0, D, numPoints,
                         m_bdata.data(),numPoints,0.0,m_dbdata.data(),numPoints);
                 }
                 break;
 
             case eModified_A:
 
                 // Note the following packing deviates from the
                 // definition in the Book by Karniadakis in that we
                 // put the vertex degrees of freedom at the lower
                 // index range to follow a more hierarchic structure.
 
                 for(i = 0; i < numPoints; ++i)
                 {
                     m_bdata[i] = 0.5*(1-z[i]);
                     m_bdata[numPoints + i] = 0.5*(1+z[i]);
                 }
 
                 mode = m_bdata.data() + 2*numPoints;
 
                 for(p = 2; p < numModes; ++p, mode += numPoints)
                 {
                     Polylib::jacobfd(numPoints, z.data(), mode, NULL, p-2,1.0,1.0);
 
                     for(i = 0; i < numPoints; ++i)
                     {
                         mode[i] *= m_bdata[i]*m_bdata[numPoints+i];
                     }
                 }
 
                 // define derivative basis
                 Blas::Dgemm('n','n',numPoints,numModes,numPoints,1.0,D,
                     numPoints,m_bdata.data(),numPoints,0.0,m_dbdata.data(),
                     numPoints);
                 break;
 
             case eModified_B:
                 {
 
                     // Note the following packing deviates from the
                     // definition in the Book by Karniadakis in two
                     // ways. 1) We put the vertex degrees of freedom
                     // at the lower index range to follow a more
                     // hierarchic structure. 2) We do not duplicate
                     // the singular vertex definition so that only a
                     // triangular number (i.e. (modes)*(modes+1)/2) of
                     // modes are required consistent with the
                     // orthogonal basis.
 
                     // In the current structure the q index runs
                     // faster than the p index so that the matrix has
                     // a more compact structure
 
                     const NekDouble *one_m_z_pow, *one_p_z;
 
                     // bdata should be of size order*(order+1)/2*zorder
 
                     // first fow
                     for(i = 0; i < numPoints; ++i)
                     {
                         m_bdata[0*numPoints + i] = 0.5*(1-z[i]);
                         m_bdata[1*numPoints + i] = 0.5*(1+z[i]);
                     }
 
                     mode = m_bdata.data() + 2*numPoints;
 
                     for(q = 2; q < numModes; ++q, mode+=numPoints)
                     {
                         Polylib::jacobfd(numPoints, z.data(), mode, NULL, q-2,1.0,1.0);
 
                         for(i = 0; i < numPoints; ++i)
                         {
                             mode[i] *= m_bdata[i]*m_bdata[numPoints+i];
                         }
                     }
 
                     // second row
                     for(i = 0; i < numPoints; ++i)
                     {
                         mode[i] = 0.5*(1-z[i]);
                     }
 
                     mode += numPoints;
 
                     for(q = 2; q < numModes; ++q, mode+=numPoints)
                     {
                         Polylib::jacobfd(numPoints, z.data(), mode, NULL, q-2,1.0,1.0);
 
                         for(i = 0; i < numPoints; ++i)
                         {
                             mode[i] *= m_bdata[i]*m_bdata[numPoints+i];
                         }
                     }
 
                     // third and higher rows
                     one_m_z_pow = m_bdata.data();
                     one_p_z     = m_bdata.data()+numPoints;
 
                     for(p = 2; p < numModes; ++p)
                     {
                         for(i = 0; i < numPoints; ++i)
                         {
                             mode[i] = m_bdata[i]*one_m_z_pow[i];
                         }
 
                         one_m_z_pow  = mode;
                         mode        += numPoints;
 
                         for(q = 1; q < numModes-p; ++q, mode+=numPoints)
                         {
                             Polylib::jacobfd(numPoints,z.data(),mode,NULL,q-1,2*p-1,1.0);
 
                             for(i = 0; i <  numPoints; ++i)
                             {
                                 mode[i] *= one_m_z_pow[i]*one_p_z[i];
                             }
                         }
                     }
 
                     Blas::Dgemm('n','n',numPoints,numModes*(numModes+1)/2,
                         numPoints,1.0,D,numPoints,
                         m_bdata.data(),numPoints,0.0,m_dbdata.data(),numPoints);
                 }
                 break;
 
 
             case eModified_C:
                 {
                     // Note the following packing deviates from the
                     // definition in the Book by Karniadakis &
                     // Sherwin in two ways. 1) We put the vertex
                     // degrees of freedom at the lower index range to
                     // follow a more hierarchic structure. 2) We do
                     // not duplicate the singular vertex definition
                     // (or the duplicated face information in the book
                     // ) so that only a tetrahedral number
                     // (i.e. (modes)*(modes+1)*(modes+2)/6) of modes
                     // are required consistent with the orthogonal
                     // basis.
 
                     // In the current structure the r index runs
                     // fastest rollowed by q and than the p index so
                     // that the matrix has a more compact structure
 
                     // Note that eModified_C is a re-organisation/
                     // duplication of eModified_B so will get a
                     // temporary Modified_B expansion and copy the
                     // correct components.
 
                     // Generate Modified_B basis;
                     BasisKey ModBKey(eModified_B,m_basisKey.GetNumModes(),
                                     m_basisKey.GetPointsKey());
                     BasisSharedPtr  ModB = BasisManager()[ModBKey];
 
                     Array<OneD, const NekDouble> ModB_data = ModB->GetBdata();
 
                     // Copy Modified_B basis into first
                     // (numModes*(numModes+1)/2)*numPoints entires of
                     // bdata.  This fills in the complete (r,p) face.
 
                     // Set up \phi^c_{p,q,r} = \phi^b_{p+q,r}
 
                     int N;
                     int B_offset = 0;
                     int offset = 0;
                     for(p = 0; p < numModes; ++p)
                     {
                         N = numPoints*(numModes-p)*(numModes-p+1)/2;
                         Vmath::Vcopy(N, &ModB_data[0]+B_offset,1,&m_bdata[0] + offset,1);
                         B_offset += numPoints*(numModes-p);
                         offset   += N;
                     }
 
                     // set up derivative of basis.
                     Blas::Dgemm('n','n',numPoints,
                                 numModes*(numModes+1)*(numModes+2)/6,
                                 numPoints,1.0,D,numPoints,
                                 m_bdata.data(),numPoints,0.0,
                                 m_dbdata.data(),numPoints);
                 }
                 break;
 
             case eModifiedPyr_C:
                 {
                     // Note the following packing deviates from the
                     // definition in the Book by Karniadakis &
                     // Sherwin in two ways. 1) We put the vertex
                     // degrees of freedom at the lower index range to
                     // follow a more hierarchic structure. 2) We do
                     // not duplicate the singular vertex definition
                     //  so that only a pyramidic  number
                     // (i.e. (modes)*(modes+1)*(2*modes+1)/6) of modes
                     // are required consistent with the orthogonal
                     // basis.
 
                     // In the current structure the r index runs
                     // fastest rollowed by q and than the p index so
                     // that the matrix has a more compact structure
 
                     // Generate Modified_B basis;
                     BasisKey ModBKey(eModified_B,m_basisKey.GetNumModes(),
                                     m_basisKey.GetPointsKey());
                     BasisSharedPtr  ModB = BasisManager()[ModBKey];
 
                     Array<OneD, const NekDouble> ModB_data = ModB->GetBdata();
 
                     // Copy Modified_B basis into first
                     // (numModes*(numModes+1)/2)*numPoints entires of
                     // bdata.
 
 
                     int N;
                     int B_offset = 0;
                     int offset = 0;
 
                     // Vertex 0,3,4, edges 3,4,7, face 4
                     N = numPoints*(numModes)*(numModes+1)/2;
                     Vmath::Vcopy(N, &ModB_data[0],1,&m_bdata[0],1);
                     offset   += N;
 
                     B_offset += numPoints*(numModes);
                     // Vertex 1 edges 5
                     N = numPoints*(numModes-1);
                     Vmath::Vcopy(N, &ModB_data[0]+B_offset,1,&m_bdata[0]+offset,1);
                     offset   += N;
 
                     // Vertex 2 edges 1,6, face 2
                     N = numPoints*(numModes-1)*(numModes)/2;
                     Vmath::Vcopy(N, &ModB_data[0]+B_offset,1,&m_bdata[0]+offset,1);
                     offset   += N;
                                  
                     B_offset += numPoints*(numModes-1);
 
                     NekDouble *one_m_z_pow, *one_p_z;
                     NekDouble *mode; 
 
                     mode = m_bdata.data() + offset;
 
                     for(p = 2; p < numModes; ++p)
                     {
                         // edges 0 2, faces 1 3
                         N = numPoints*(numModes-p);
                         Vmath::Vcopy(N, &ModB_data[0]+B_offset,1,mode,1);
                         mode   += N;
                         Vmath::Vcopy(N, &ModB_data[0]+B_offset,1,mode,1);
                         mode   += N;
                         B_offset += N;
                     
                         one_p_z     = m_bdata.data()+numPoints;
 
                         for(q = 2; q < numModes; ++q)
                         {
                             // face 0 
                             for(i = 0; i < numPoints; ++i)
                             {
                                 // [(1-z)/2]^{p+q-2} Note in book it
                                 // seems to suggest p+q-1 but that
                                 // does not seem to give complete
                                 // polynomial space for pyramids
                                 mode[i] = pow(m_bdata[i],p+q-2);
                             }
 
                             one_m_z_pow  = mode;
                             mode        += numPoints;
                             
                             // interior 
                             for(int r = 1; r < numModes-std::max(p,q); ++r)
                             {
                                 Polylib::jacobfd(numPoints,z.data(),mode,NULL,r-1,2*p+2*q-3,1.0);
                                 
                                 for(i = 0; i <  numPoints; ++i)
                                 {
                                     mode[i] *= one_m_z_pow[i]*one_p_z[i];
                                 }
                                 mode   += numPoints;
                             }
                         }
                     
                     }
 
                     // set up derivative of basis.
                     Blas::Dgemm('n','n',numPoints,
                                 numModes*(numModes+1)*(2*numModes+1)/6,
                                 numPoints,1.0,D,numPoints,
                                 m_bdata.data(),numPoints,0.0,
                                 m_dbdata.data(),numPoints);
                 }
                 break;
 
             case eGLL_Lagrange:
                 {
                     mode = m_bdata.data();
                     std::shared_ptr< Points<NekDouble> > m_points = PointsManager()[PointsKey(numModes, eGaussLobattoLegendre)];
                     const Array<OneD, const NekDouble>& zp(m_points->GetZ());
 
                     for (p=0; p<numModes; ++p, mode += numPoints)
                     {
                         for(q = 0; q < numPoints; ++q)
                         {
                             mode[q] = Polylib::hglj(p, z[q], zp.data(), numModes, 0.0, 0.0);
                         }
                     }
 
                     // define derivative basis
                     Blas::Dgemm('n', 'n', numPoints, numModes, numPoints, 1.0,
                                 D, numPoints, m_bdata.data(), numPoints, 0.0,
                                 m_dbdata.data(), numPoints);
                     
                 }//end scope
                 break;
             case eGauss_Lagrange:
                 {
                     mode = m_bdata.data();
                     std::shared_ptr< Points<NekDouble> > m_points = PointsManager()[PointsKey(numModes, eGaussGaussLegendre)];
                     const Array<OneD, const NekDouble>& zp(m_points->GetZ());
                     
                     for (p=0; p<numModes; ++p,mode += numPoints)
                     {
                         for(q = 0; q < numPoints; ++q)
                         {
                             mode[q] = Polylib::hgj(p, z[q], zp.data(), numModes, 0.0, 0.0);
                         }
                     }
                     
                     // define derivative basis
                     Blas::Dgemm('n', 'n', numPoints, numModes, numPoints, 1.0,
                                 D, numPoints, m_bdata.data(), numPoints, 0.0,
                                 m_dbdata.data(), numPoints);
                     
                 }//end scope
                 break;
             case eFourier:
 
                 ASSERTL0(numModes%2==0, "Fourier modes should be a factor of 2");
 
                 for(i = 0; i < numPoints; ++i)
                 {
                     m_bdata[i] = 1.0;
                     m_bdata[numPoints+i] = 0.0;
 
                     m_dbdata[i] = m_dbdata[numPoints+i] = 0.0;
                 }
 
                 for (p=1; p < numModes/2; ++p)
                 {
                     for(i = 0; i < numPoints; ++i)
                     {
                         m_bdata[ 2*p   *numPoints+i] =  cos(p*M_PI* (z[i]+1) );
                         m_bdata[(2*p+1)*numPoints+i] = -sin(p*M_PI* (z[i]+1) );
 
                         m_dbdata[ 2*p   *numPoints+i] = -p*M_PI*sin(p*M_PI* (z[i]+1) );
                         m_dbdata[(2*p+1)*numPoints+i] = -p*M_PI*cos(p*M_PI* (z[i]+1) );
                     }
                 }
 
                 break;
                     
             
             // Fourier Single Mode (1st mode)
             case eFourierSingleMode:
                     
                 for(i = 0; i < numPoints; ++i)
                 {
                     m_bdata[i] = cos(M_PI* (z[i]+1) );
                     m_bdata[numPoints+i] = -sin(M_PI* (z[i]+1) );
                         
                     m_dbdata[i] = -M_PI*sin(M_PI* (z[i]+1) );
                     m_dbdata[numPoints+i] = -M_PI*cos(M_PI* (z[i]+1) );
                 }
                     
                 for (p=1; p < numModes/2; ++p)
                 {
                     for(i = 0; i < numPoints; ++i)
                     {
                         m_bdata[ 2*p   *numPoints+i] = 0.;
                         m_bdata[(2*p+1)*numPoints+i] = 0.;
                             
                         m_dbdata[ 2*p   *numPoints+i] = 0.;
                         m_dbdata[(2*p+1)*numPoints+i] = 0.;
                     }
                 }
                 break;
 
             //Fourier Real Half Mode
             case eFourierHalfModeRe:
                 m_bdata[0] = cos(M_PI*z[0]);
                 m_dbdata[0] = -M_PI*sin(M_PI*z[0]);
                 break;
 
             //Fourier Imaginary Half Mode
             case eFourierHalfModeIm:
                 m_bdata[0] = -sin(M_PI*z[0]);
                 m_dbdata[0] = -M_PI*cos(M_PI*z[0]);
                 break;
 
             case eChebyshev:
                 {
                     mode = m_bdata.data();
 
                     for (p=0,scal = 1; p<numModes; ++p,mode += numPoints)
                     {
                         Polylib::jacobfd(numPoints, z.data(), mode, NULL, p, -0.5, -0.5);
 
                         for(i = 0; i < numPoints; ++i)
                         {
                             mode[i] *= scal;
                         }
 
                         scal *= 4*(p+1)*(p+1)/(2*p+2)/(2*p+1);
                     }
 
                     // Define derivative basis
                     Blas::Dgemm('n', 'n', numPoints, numModes, numPoints, 1.0,
                                 D, numPoints, m_bdata.data(), numPoints, 0.0,
                                 m_dbdata.data(), numPoints);
                 }
                 break;
 
             case eMonomial:
                 {
                     int P = numModes - 1;
                     NekDouble *mode = m_bdata.data();
 
                     for( int p = 0; p <= P; ++p, mode += numPoints )
                     {
                         for( int i = 0; i < numPoints; ++i )
                         {
                             mode[i] = pow(z[i], p);
                         }
                     }
 
                     // define derivative basis
                     Blas::Dgemm('n', 'n', numPoints, numModes, numPoints, 1.0,
                                 D, numPoints, m_bdata.data(), numPoints, 0.0, 
                                 m_dbdata.data(),numPoints);
                 }//end scope
                 break;
             default:
                 NEKERROR(ErrorUtil::efatal, "Basis Type not known or "
                                 "not implemented at this time.");
             }
         }

◆ GetBasisKey()

const BasisKey Nektar::LibUtilities::Basis::GetBasisKey ( ) const

inline

Returns the specification for the Basis.

Definition at line 325 of file Basis.h.

References LIB_UTILITIES_EXPORT.

Referenced by export_Basis().

             {
                 return m_basisKey;
             }

◆ GetBasisType()

BasisType Nektar::LibUtilities::Basis::GetBasisType ( ) const

inline

Return the type of expansion basis.

Definition at line 247 of file Basis.h.

Referenced by export_Basis(), and GenBasis().

             {
                 return m_basisKey.GetBasisType();
             }

◆ GetBdata()

const Array<OneD, const NekDouble>& Nektar::LibUtilities::Basis::GetBdata ( ) const

inline

Return basis definition array m_bdata.

Definition at line 313 of file Basis.h.

Referenced by Nektar::LocalRegions::Expansion1D::AddHDGHelmholtzTraceTerms(), Nektar::LocalRegions::Expansion1D::AddNormTraceInt(), and export_Basis().

             {
                 return m_bdata;
             }

◆ GetD()

const std::shared_ptr<NekMatrix<NekDouble> >& Nektar::LibUtilities::Basis::GetD ( Direction dir = xDir ) const

inline

Definition at line 280 of file Basis.h.

             {
                 return m_points->GetD(dir);
             }

◆ GetDbdata()

const Array<OneD, const NekDouble>& Nektar::LibUtilities::Basis::GetDbdata ( ) const

inline

Return basis definition array m_dbdata.

Definition at line 319 of file Basis.h.

Referenced by export_Basis().

             {
                 return m_dbdata;
             }

◆ GetI() [1/2]

const std::shared_ptr<NekMatrix<NekDouble> > Nektar::LibUtilities::Basis::GetI ( const Array< OneD, const NekDouble > & x )

inline

Definition at line 286 of file Basis.h.

             {
                 return m_points->GetI(x);
             }

◆ GetI() [2/2]

const std::shared_ptr<NekMatrix<NekDouble> > Nektar::LibUtilities::Basis::GetI ( const BasisKey & bkey )

inline

Definition at line 292 of file Basis.h.

References ASSERTL0, Nektar::LibUtilities::PointsKey::GetPointsDim(), and Nektar::LibUtilities::BasisKey::GetPointsKey().

             {
                 ASSERTL0(bkey.GetPointsKey().GetPointsDim()==1,
                          "Interpolation only to other 1d basis");
                 return m_InterpManager[bkey];
             }

◆ GetNumModes()

int Nektar::LibUtilities::Basis::GetNumModes ( ) const

inline

Return order of basis from the basis specification.

Definition at line 223 of file Basis.h.

Referenced by Nektar::LibUtilities::BasisKey::Collocation(), Nektar::LibUtilities::BasisKey::ExactIprodInt(), export_Basis(), GenBasis(), Initialize(), and main().

             {
                 return m_basisKey.GetNumModes();
             }

◆ GetNumPoints()

int Nektar::LibUtilities::Basis::GetNumPoints ( ) const

inline

Return the number of points from the basis specification.

Definition at line 235 of file Basis.h.

Referenced by Nektar::LibUtilities::BasisKey::Collocation(), Nektar::LibUtilities::BasisKey::ExactIprodInt(), export_Basis(), and GenBasis().

             {
                 return m_basisKey.GetNumPoints();
             }

◆ GetPointsKey()

PointsKey Nektar::LibUtilities::Basis::GetPointsKey ( ) const

inline

Return the points specification for the basis.

Definition at line 253 of file Basis.h.

Referenced by export_Basis().

             {
                 return m_basisKey.GetPointsKey();
             }

◆ GetPointsType()

PointsType Nektar::LibUtilities::Basis::GetPointsType ( ) const

inline

Return the type of quadrature.

Definition at line 259 of file Basis.h.

Referenced by Nektar::LibUtilities::BasisKey::Collocation(), and Nektar::LibUtilities::BasisKey::ExactIprodInt().

             {
                 return m_basisKey.GetPointsType();
             }

◆ GetTotNumModes()

int Nektar::LibUtilities::Basis::GetTotNumModes ( ) const

inline

Return total number of modes from the basis specification.

Definition at line 229 of file Basis.h.

Referenced by export_Basis().

             {
                 return m_basisKey.GetTotNumModes();
             }

◆ GetTotNumPoints()

int Nektar::LibUtilities::Basis::GetTotNumPoints ( ) const

inline

Return total number of points from the basis specification.

Definition at line 241 of file Basis.h.

Referenced by export_Basis(), and Initialize().

             {
                 return m_basisKey.GetTotNumPoints();
             }

◆ GetW()

const Array<OneD, const NekDouble>& Nektar::LibUtilities::Basis::GetW ( ) const

inline

Definition at line 269 of file Basis.h.

Referenced by export_Basis().

             {
                 return m_points->GetW();
             }

◆ GetZ()

const Array<OneD, const NekDouble>& Nektar::LibUtilities::Basis::GetZ ( ) const

inline

Definition at line 264 of file Basis.h.

Referenced by export_Basis().

             {
                 return m_points->GetZ();
             }

◆ GetZW()

void Nektar::LibUtilities::Basis::GetZW	(	Array< OneD, const NekDouble > &	z,
		Array< OneD, const NekDouble > &	w
	)		const

inline

Definition at line 274 of file Basis.h.

             {
                 m_points->GetZW(z,w);
             }

◆ Initialize()

void Nektar::LibUtilities::Basis::Initialize ( )

virtual

Definition at line 112 of file Foundations/Basis.cpp.

References ASSERTL0, GenBasis(), GetNumModes(), and GetTotNumPoints().

Referenced by export_Basis().

         {
             ASSERTL0(GetNumModes()>0, "Cannot call Basis initialisation with zero or negative order");
             ASSERTL0(GetTotNumPoints()>0, "Cannot call Basis initialisation with zero or negative numbers of points");
 
             GenBasis();
         }

Member Data Documentation

◆ initBasisManager

bool Nektar::LibUtilities::Basis::initBasisManager

staticprivate

Initial value:

= {
            BasisManager().RegisterGlobalCreator(Basis::Create)
        }

Definition at line 341 of file Basis.h.

◆ m_basisKey

BasisKey Nektar::LibUtilities::Basis::m_basisKey

protected

Basis specification.

Definition at line 333 of file Basis.h.

Referenced by CalculateInterpMatrix(), and GenBasis().

◆ m_bdata

Array<OneD, NekDouble> Nektar::LibUtilities::Basis::m_bdata

protected

Basis definition.

Definition at line 335 of file Basis.h.

Referenced by GenBasis().

◆ m_dbdata

Array<OneD, NekDouble> Nektar::LibUtilities::Basis::m_dbdata

protected

Derivative Basis definition.

Definition at line 336 of file Basis.h.

Referenced by GenBasis().

◆ m_InterpManager

NekManager<BasisKey, NekMatrix<NekDouble>, BasisKey::opLess> Nektar::LibUtilities::Basis::m_InterpManager

protected

Definition at line 338 of file Basis.h.

Referenced by Basis().

◆ m_points

PointsSharedPtr Nektar::LibUtilities::Basis::m_points

protected

Set of points.

Definition at line 334 of file Basis.h.

Referenced by GenBasis().

Public Member Functions

Static Public Member Functions

Protected Attributes

Private Member Functions

Static Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ ~Basis()

◆ Basis() [1/2]

◆ Basis() [2/2]

Member Function Documentation

◆ CalculateInterpMatrix()

◆ Collocation()

◆ Create()

◆ ExactIprodInt()

◆ GenBasis()

◆ GetBasisKey()

◆ GetBasisType()

◆ GetBdata()

◆ GetD()

◆ GetDbdata()

◆ GetI() [1/2]

◆ GetI() [2/2]

◆ GetNumModes()

◆ GetNumPoints()

◆ GetPointsKey()

◆ GetPointsType()

◆ GetTotNumModes()

◆ GetTotNumPoints()

◆ GetW()

◆ GetZ()

◆ GetZW()

◆ Initialize()

Member Data Documentation

◆ initBasisManager

◆ m_basisKey

◆ m_bdata

◆ m_dbdata

◆ m_InterpManager

◆ m_points