Nektar::LibUtilities::Basis Class Reference

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

#include <Basis.h>

Public Member Functions
virtual	~Basis ()
	Destructor.
int	GetNumModes () const
	Return order of basis from the basis specification.
int	GetTotNumModes () const
	Return total number of modes from the basis specification.
int	GetNumPoints () const
	Return the number of points from the basis specification.
int	GetTotNumPoints () const
	Return total number of points from the basis specification.
BasisType	GetBasisType () const
	Return the type of expansion basis.
PointsKey	GetPointsKey () const
	Return the points specification for the basis.
PointsType	GetPointsType () const
	Return the type of quadrature.
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.
bool	Collocation () const
	Determine if basis has collocation properties.
const Array< OneD, const NekDouble > &	GetBdata () const
	Return basis definition array m_bdata.
const Array< OneD, const NekDouble > &	GetDbdata () const
	Return basis definition array m_dbdata.
const BasisKey	GetBasisKey () const
	Returns the specification for the Basis.
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.

Protected Attributes
BasisKey	m_basisKey
	Basis specification.
PointsSharedPtr	m_points
	Set of points.
Array< OneD, NekDouble >	m_bdata
	Basis definition.
Array< OneD, NekDouble >	m_dbdata
	Derivative Basis definition.
NekManager< BasisKey, NekMatrix< NekDouble >, BasisKey::opLess >	m_InterpManager

Private Member Functions
	Basis (const BasisKey &bkey)
	Private constructor with BasisKey.
	Basis ()=delete
	Private default constructor.
std::shared_ptr< NekMatrix< NekDouble > >	CalculateInterpMatrix (const BasisKey &tbasis0)
	Calculate the interpolation Matrix for coefficient from one base (m_basisKey) to another (tbasis0) In the new version, we do not enforce the interpolation matrix to be square. The output dimension will be determined by tBasis and input dimension will be determined by m_basisKey.
void	GenBasis ()
	Generate appropriate basis and their derivatives.

Static Private Attributes
static bool	initBasisManager

Detailed Description

Represents a basis of a given type.

Definition at line 197 of file Basis.h.

Constructor & Destructor Documentation

~Basis()

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

inlinevirtual

Destructor.

Definition at line 204 of file Basis.h.

{};

Basis() [1/2]

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

private

Private constructor with BasisKey.

Definition at line 90 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 ( )

privatedelete

Private default constructor.

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) In the new version, we do not enforce the interpolation matrix to be square. The output dimension will be determined by tBasis and input dimension will be determined by m_basisKey.

• If tBasis.GetNumModes() > m_basisKey.GetNumModes(), the interpolation matrix will be $B^{(l-1,l)T}_f B^{(-l)T}_t$, using tBasis.GetNumModes() as basis points.
• If tBasis.GetNumModes() <= m_basisKey.GetNumModes(), the interpolation matrix will be $B^{(-l)}_f B^{(l,l-1)}_t$, using m_basisKey.GetNumModes() as basis points. which is the same as the old version.

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

{
    // Get the number of modes of the target basis
    size_t tdim = tbasis0.GetNumModes();
    // Get the number of modes of the source basis
    size_t fdim = m_basisKey.GetNumModes();
    if (tdim < fdim)
    {
        size_t dim = fdim;
        const PointsKey pkey(dim, LibUtilities::eGaussLobattoLegendre);
        BasisKey fbkey(m_basisKey.GetBasisType(), fdim, pkey);
        BasisKey tbkey(tbasis0.GetBasisType(), tdim, 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(fdim, dim, fB_data);
        NekMatrix<NekDouble> tB(dim, tdim, tB_data);
 
        // Invert the "from" matrix
        fB.Invert();
 
        // Compute transformation matrix
        Array<OneD, NekDouble> zero1D(tdim * fdim, 0.0);
        std::shared_ptr<NekMatrix<NekDouble>> ftB(
            MemoryManager<NekMatrix<NekDouble>>::AllocateSharedPtr(fdim, tdim,
                                                                   zero1D));
        // tu := tB^T * fB^(-1)^T  fu = ftB fu
        (*ftB) = fB * tB;
        // this operation will not acturally do the transpose, just set the
        // transpose flag to 'T' ftB->Transpose(); we need to do the transpose
        // manually
        std::shared_ptr<NekMatrix<NekDouble>> ftB_T(
            MemoryManager<NekMatrix<NekDouble>>::AllocateSharedPtr(tdim, fdim,
                                                                   0.0));
        for (size_t i = 0; i < tdim; ++i)
        {
            for (size_t j = 0; j < fdim; ++j)
            {
                (*ftB_T)(i, j) = (*ftB)(j, i);
            }
        }
        return ftB_T;
    }
    else // original version
    {
        size_t dim = tdim;
        const PointsKey pkey(dim, LibUtilities::eGaussLobattoLegendre);
        BasisKey fbkey(m_basisKey.GetBasisType(), fdim, pkey);
        BasisKey tbkey(tbasis0.GetBasisType(), tdim, 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, fdim, fB_data);
        NekMatrix<NekDouble> tB(dim, tdim, tB_data);
 
        // Invert the "to" matrix: tu := tB^(-1)*fB fu = ftB fu
        tB.Invert();
 
        // Compute transformation matrix
        Array<OneD, NekDouble> zero1D(tdim * fdim, 0.0);
        std::shared_ptr<NekMatrix<NekDouble>> ftB(
            MemoryManager<NekMatrix<NekDouble>>::AllocateSharedPtr(tdim, fdim,
                                                                   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 295 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 99 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 ( ) const

inline

Determine if basis has exact integration for inner product.

Definition at line 289 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 276 of file Foundations/Basis.cpp.

{
    size_t i, p, q;
    NekDouble scal;
    Array<OneD, NekDouble> modeSharedArray;
    NekDouble *mode;
    Array<OneD, const NekDouble> z;
    Array<OneD, const NekDouble> w;
 
    m_points->GetZW(z, w);
 
    const NekDouble *D = &(m_points->GetD()->GetPtr())[0];
    size_t numModes    = GetNumModes();
    size_t 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, nullptr, 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);
        } // end scope
        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 Spencer'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:
        {
            mode = m_bdata.data();
 
            for (size_t p = 0; p < numModes; ++p)
            {
                for (size_t q = 0; q < numModes - p; ++q, mode += numPoints)
                {
                    Polylib::jacobfd(numPoints, z.data(), mode, nullptr, q,
                                     2 * p + 1.0, 0.0);
                    for (size_t 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 Spencer'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:
        {
            size_t P = numModes - 1, Q = numModes - 1, R = numModes - 1;
            mode = m_bdata.data();
 
            for (size_t p = 0; p <= P; ++p)
            {
                for (size_t q = 0; q <= Q - p; ++q)
                {
                    for (size_t r = 0; r <= R - p - q; ++r, mode += numPoints)
                    {
                        Polylib::jacobfd(numPoints, z.data(), mode, nullptr, r,
                                         2 * p + 2 * q + 2.0, 0.0);
                        for (size_t 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:
        {
            size_t P = numModes - 1, Q = numModes - 1, R = numModes - 1;
            mode = m_bdata.data();
 
            for (size_t p = 0; p <= P; ++p)
            {
                for (size_t q = 0; q <= Q; ++q)
                {
                    for (size_t 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.
                        // size_t pq = max(p + q -1,0);
                        size_t pq = std::max(p + q, size_t(0));
 
                        Polylib::jacobfd(numPoints, z.data(), mode, nullptr, r,
                                         2 * pq + 2.0, 0.0);
                        for (size_t 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, nullptr, 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, nullptr, 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, nullptr, 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, nullptr, 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}
 
            size_t N;
            size_t B_offset = 0;
            size_t 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.
 
            size_t N;
            size_t B_offset = 0;
            size_t 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;
 
            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 (size_t r = 1; r < numModes - std::max(p, q); ++r)
                    {
                        Polylib::jacobfd(numPoints, z.data(), mode, nullptr,
                                         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);
        } // end scope
        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));
                }
            }
        } // end scope
        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.;
                }
            }
        } // end scope
        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]);
        } // end scope
        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]);
        } // end scope
        break;
 
        case eChebyshev:
        {
            mode = m_bdata.data();
 
            for (p = 0, scal = 1; p < numModes; ++p, mode += numPoints)
            {
                Polylib::jacobfd(numPoints, z.data(), mode, nullptr, 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);
        } // end scope
        break;
 
        case eMonomial:
        {
            mode = m_bdata.data();
 
            for (size_t p = 0; p < numModes; ++p, mode += numPoints)
            {
                for (size_t 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::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, 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 264 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 313 of file Basis.h.

    {
        return m_basisKey;
    }

References m_basisKey.

Referenced by export_Basis().

GetBasisType()

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

inline

Return the type of expansion basis.

Definition at line 231 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 301 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 269 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 307 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 275 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 281 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 207 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 219 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 237 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 243 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 213 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 225 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 253 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 248 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 258 of file Basis.h.

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

References m_points.

Initialize()

void Nektar::LibUtilities::Basis::Initialize

(

)

Definition at line 107 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 329 of file Basis.h.

m_basisKey

BasisKey Nektar::LibUtilities::Basis::m_basisKey

protected

Basis specification.

Definition at line 321 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 323 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 324 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 326 of file Basis.h.

Referenced by Basis(), and GetI().

m_points

PointsSharedPtr Nektar::LibUtilities::Basis::m_points

protected

Set of points.

Definition at line 322 of file Basis.h.

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