Foundations/Basis.cpp

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File: Basis.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Basis definition
//
///////////////////////////////////////////////////////////////////////////////
 
#include <LibUtilities/Foundations/Basis.h>
#include <LibUtilities/Foundations/ManagerAccess.h>
#include <LibUtilities/Polylib/Polylib.h>
 
namespace Nektar::LibUtilities
{
bool Basis::initBasisManager = {
    BasisManager().RegisterGlobalCreator(Basis::Create)};
 
bool operator<(const BasisKey &lhs, const BasisKey &rhs)
{
    PointsKey lhsPointsKey = lhs.GetPointsKey();
    PointsKey rhsPointsKey = rhs.GetPointsKey();
 
    if (lhsPointsKey < rhsPointsKey)
    {
        return true;
    }
    if (lhsPointsKey != rhsPointsKey)
    {
        return false;
    }
 
    if (lhs.m_nummodes < rhs.m_nummodes)
    {
        return true;
    }
    if (lhs.m_nummodes > rhs.m_nummodes)
    {
        return false;
    }
 
    return (lhs.m_basistype < rhs.m_basistype);
}
 
bool operator>(const BasisKey &lhs, const BasisKey &rhs)
{
    return (rhs < lhs);
}
 
bool BasisKey::opLess::operator()(const BasisKey &lhs,
                                  const BasisKey &rhs) const
{
    return (lhs.m_basistype < rhs.m_basistype);
}
 
std::ostream &operator<<(std::ostream &os, const BasisKey &rhs)
{
    os << "NumModes: " << rhs.GetNumModes()
       << " BasisType: " << BasisTypeMap[rhs.GetBasisType()];
    os << " " << rhs.GetPointsKey() << std::endl;
 
    return os;
}
 
Basis::Basis(const BasisKey &bkey)
    : 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));
}
 
std::shared_ptr<Basis> Basis::Create(const BasisKey &bkey)
{
    std::shared_ptr<Basis> returnval(new Basis(bkey));
    returnval->Initialize();
 
    return returnval;
}
 
void Basis::Initialize()
{
    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();
}
 
/** \brief 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.
 */
std::shared_ptr<NekMatrix<NekDouble>> Basis::CalculateInterpMatrix(
    const BasisKey &tbasis0)
{
    // 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;
    }
}
 
// Method used to generate appropriate basis
/** The following expansions are generated depending on the
 * enum type defined in \a 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
 *
 * \a eModified_A:
 *
 * m_bdata[i + j*m_numpoints] =
 * \f$ \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 . \f$
 *
 * \a eModified_B:
 *
 * m_bdata[n(i,j) + k*m_numpoints] =
 * \f$ \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 . , \f$
 *
 * where \f$ n(i,j) \f$ is a consecutive ordering of the
 * triangular indices \f$ 0 \leq i, i+j < P \f$ where \a j
 * runs fastest.
 *
 *
 * \a eModified_C:
 *
 * m_bdata[n(i,j,k) + l*m_numpoints] =
 * \f$ \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 . , \f$
 *
 * where \f$ n(i,j,k) \f$ is a consecutive ordering of the
 * triangular indices \f$ 0 \leq i, i+j, i+j+k < P \f$ where \a k
 * runs fastest, then \a j and finally \a i.
 *
 */
void Basis::GenBasis()
{
    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.");
    }
}
 
/** \brief Determine if polynomial basis can be exactly integrated
 *  with itself
 */
bool BasisKey::ExactIprodInt(void) const
{
    bool returnval = false;
 
    switch (GetPointsType())
    {
        case eGaussGaussLegendre:
        case eGaussRadauMLegendre:
        case eGaussRadauPLegendre:
        case eGaussLobattoLegendre:
        case eGaussKronrodLegendre:
        case eGaussRadauKronrodMLegendre:
        case eGaussLobattoKronrodLegendre:
            returnval = (GetNumPoints() >= GetNumModes());
            break;
 
        default:
            break;
    }
 
    return returnval;
}
 
/** \brief Determine if basis has collocation property,
 *  i.e. GLL_Lagrange with Lobatto integration of appropriate order,
 *  Gauss_Lagrange with Gauss integration of appropriate order.
 */
bool BasisKey::Collocation() const
{
    return ((m_basistype == eGLL_Lagrange &&
             GetPointsType() == eGaussLobattoLegendre &&
             GetNumModes() == GetNumPoints()) ||
            (m_basistype == eGauss_Lagrange &&
             GetPointsType() == eGaussGaussLegendre &&
             GetNumModes() == GetNumPoints()));
}
 
// BasisKey compared to BasisKey
bool operator==(const BasisKey &x, const BasisKey &y)
{
    return (x.GetPointsKey() == y.GetPointsKey() &&
            x.m_basistype == y.m_basistype &&
            x.GetNumModes() == y.GetNumModes());
}
 
// BasisKey* compared to BasisKey
bool operator==(const BasisKey *x, const BasisKey &y)
{
    return (*x == y);
}
 
// \brief BasisKey compared to BasisKey*
bool operator==(const BasisKey &x, const BasisKey *y)
{
    return (x == *y);
}
 
// \brief BasisKey compared to BasisKey
bool operator!=(const BasisKey &x, const BasisKey &y)
{
    return (!(x == y));
}
 
// BasisKey* compared to BasisKey
bool operator!=(const BasisKey *x, const BasisKey &y)
{
    return (!(*x == y));
}
 
// BasisKey compared to BasisKey*
bool operator!=(const BasisKey &x, const BasisKey *y)
{
    return (!(x == *y));
}
 
} // namespace Nektar::LibUtilities