Nektar::LibUtilities::Basis Class Reference

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 Array< OneD, const NekDouble > &	GetBaryWeights () 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...
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 212 of file Basis.h.

Constructor & Destructor Documentation

~Basis()

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

inlinevirtual

Destructor.

Definition at line 219 of file Basis.h.

{};

Basis() [1/2]

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

private

Private constructor with BasisKey.

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

    : 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));
}

References CalculateInterpMatrix(), and m_InterpManager.

Basis() [2/2]

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

inlineprivate

Private default constructor.

Definition at line 350 of file Basis.h.

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

References Nektar::ErrorUtil::efatal, and NEKERROR.

Referenced by Create().

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 127 of file Foundations/Basis.cpp.

{
    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;
}

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

Referenced by Basis().

Collocation()

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

inline

Determine if basis has collocation properties.

Definition at line 310 of file Basis.h.

    {
        return m_basisKey.Collocation();
    }

References Nektar::LibUtilities::BasisKey::Collocation(), and m_basisKey.

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 105 of file Foundations/Basis.cpp.

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

References Basis().

ExactIprodInt()

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

inline

Determine if basis has exact integration for inner product.

Definition at line 304 of file Basis.h.

    {
        return m_basisKey.ExactIprodInt();
    }

References Nektar::LibUtilities::BasisKey::ExactIprodInt(), and m_basisKey.

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 223 of file Foundations/Basis.cpp.

{
    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.");
    }
}

References ASSERTL0, Nektar::LibUtilities::BasisManager(), Blas::Dgemm(), Nektar::LibUtilities::eChebyshev, Nektar::ErrorUtil::efatal, Nektar::LibUtilities::eFourier, Nektar::LibUtilities::eFourierHalfModeIm, Nektar::LibUtilities::eFourierHalfModeRe, Nektar::LibUtilities::eFourierSingleMode, Nektar::LibUtilities::eGauss_Lagrange, Nektar::LibUtilities::eGaussGaussLegendre, Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eGLL_Lagrange, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::LibUtilities::eModified_C, Nektar::LibUtilities::eModifiedPyr_C, Nektar::LibUtilities::eOrtho_A, Nektar::LibUtilities::eOrtho_B, Nektar::LibUtilities::eOrtho_C, Nektar::LibUtilities::eOrthoPyr_C, GetBasisType(), Nektar::LibUtilities::BasisKey::GetNumModes(), GetNumModes(), GetNumPoints(), Nektar::LibUtilities::BasisKey::GetPointsKey(), Polylib::hgj(), Polylib::hglj(), Polylib::jacobfd(), m_basisKey, m_bdata, m_dbdata, m_points, NEKERROR, Nektar::LibUtilities::P, CellMLToNektar.cellml_metadata::p, Nektar::LibUtilities::PointsManager(), tinysimd::sqrt(), and Vmath::Vcopy().

Referenced by Initialize().

GetBaryWeights()

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

inline

Definition at line 279 of file Basis.h.

    {
        return m_points->GetBaryWeights();
    }

References m_points.

GetBasisKey()

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

inline

Returns the specification for the Basis.

Definition at line 328 of file Basis.h.

    {
        return m_basisKey;
    }

References m_basisKey.

Referenced by export_Basis(), and Nektar::MultiRegions::ExpListHomogeneous2D::GenHomogeneous2DBlockMatrix().

GetBasisType()

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

inline

Return the type of expansion basis.

Definition at line 246 of file Basis.h.

    {
        return m_basisKey.GetBasisType();
    }

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

Referenced by export_Basis(), and GenBasis().

GetBdata()

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

inline

Return basis definition array m_bdata.

Definition at line 316 of file Basis.h.

    {
        return m_bdata;
    }

References m_bdata.

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

GetD()

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

inline

Definition at line 284 of file Basis.h.

    {
        return m_points->GetD(dir);
    }

References m_points.

GetDbdata()

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

inline

Return basis definition array m_dbdata.

Definition at line 322 of file Basis.h.

    {
        return m_dbdata;
    }

References m_dbdata.

Referenced by export_Basis().

GetI() [1/2]

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

inline

Definition at line 290 of file Basis.h.

    {
        return m_points->GetI(x);
    }

References m_points.

GetI() [2/2]

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

inline

Definition at line 296 of file Basis.h.

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

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

GetNumModes()

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

inline

Return order of basis from the basis specification.

Definition at line 222 of file Basis.h.

    {
        return m_basisKey.GetNumModes();
    }

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

Referenced by export_Basis(), GenBasis(), Initialize(), and main().

GetNumPoints()

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

inline

Return the number of points from the basis specification.

Definition at line 234 of file Basis.h.

    {
        return m_basisKey.GetNumPoints();
    }

References Nektar::LibUtilities::BasisKey::GetNumPoints(), and m_basisKey.

Referenced by export_Basis(), and GenBasis().

GetPointsKey()

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

inline

Return the points specification for the basis.

Definition at line 252 of file Basis.h.

    {
        return m_basisKey.GetPointsKey();
    }

References Nektar::LibUtilities::BasisKey::GetPointsKey(), and m_basisKey.

Referenced by export_Basis().

GetPointsType()

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

inline

Return the type of quadrature.

Definition at line 258 of file Basis.h.

    {
        return m_basisKey.GetPointsType();
    }

References Nektar::LibUtilities::BasisKey::GetPointsType(), and m_basisKey.

GetTotNumModes()

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

inline

Return total number of modes from the basis specification.

Definition at line 228 of file Basis.h.

    {
        return m_basisKey.GetTotNumModes();
    }

References Nektar::LibUtilities::BasisKey::GetTotNumModes(), and m_basisKey.

Referenced by export_Basis().

GetTotNumPoints()

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

inline

Return total number of points from the basis specification.

Definition at line 240 of file Basis.h.

    {
        return m_basisKey.GetTotNumPoints();
    }

References Nektar::LibUtilities::BasisKey::GetTotNumPoints(), and m_basisKey.

Referenced by export_Basis(), and Initialize().

GetW()

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

inline

Definition at line 268 of file Basis.h.

    {
        return m_points->GetW();
    }

References m_points.

Referenced by export_Basis().

GetZ()

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

inline

Definition at line 263 of file Basis.h.

    {
        return m_points->GetZ();
    }

References m_points.

Referenced by export_Basis().

GetZW()

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

inline

Definition at line 273 of file Basis.h.

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

References m_points.

Initialize()

void Nektar::LibUtilities::Basis::Initialize

(

)

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

{
    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();
}

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

Referenced by export_Basis().

Member Data Documentation

initBasisManager

bool Nektar::LibUtilities::Basis::initBasisManager

staticprivate

Initial value:

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

Definition at line 344 of file Basis.h.

m_basisKey

BasisKey Nektar::LibUtilities::Basis::m_basisKey

protected

Basis specification.

Definition at line 336 of file Basis.h.

Referenced by CalculateInterpMatrix(), Collocation(), ExactIprodInt(), GenBasis(), GetBasisKey(), GetBasisType(), GetNumModes(), GetNumPoints(), GetPointsKey(), GetPointsType(), GetTotNumModes(), and GetTotNumPoints().

m_bdata

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

protected

Basis definition.

Definition at line 338 of file Basis.h.

Referenced by GenBasis(), and GetBdata().

m_dbdata

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

protected

Derivative Basis definition.

Definition at line 339 of file Basis.h.

Referenced by GenBasis(), and GetDbdata().

m_InterpManager

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

protected

Definition at line 341 of file Basis.h.

Referenced by Basis(), and GetI().

m_points

PointsSharedPtr Nektar::LibUtilities::Basis::m_points

protected

Set of points.

Definition at line 337 of file Basis.h.

Referenced by GenBasis(), GetBaryWeights(), GetD(), GetI(), GetW(), GetZ(), and GetZW().