StdQuadExp.cpp

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File: StdQuadExp.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: Quadrilateral routines built upon StdExpansion2D
//
///////////////////////////////////////////////////////////////////////////////
 
#include <LibUtilities/Foundations/InterpCoeff.h>
#include <LibUtilities/Foundations/ManagerAccess.h>
#include <StdRegions/StdQuadExp.h>
 
using namespace std;
 
namespace Nektar::StdRegions
{
 
/** \brief Constructor using BasisKey class for quadrature
 *  points and order definition
 */
StdQuadExp::StdQuadExp(const LibUtilities::BasisKey &Ba,
                       const LibUtilities::BasisKey &Bb)
    : StdExpansion(Ba.GetNumModes() * Bb.GetNumModes(), 2, Ba, Bb),
      StdExpansion2D(Ba.GetNumModes() * Bb.GetNumModes(), Ba, Bb)
{
}
 
/////////////////////////
// Integration Methods //
/////////////////////////
 
NekDouble StdQuadExp::v_Integral(const Array<OneD, const NekDouble> &inarray)
{
    Array<OneD, const NekDouble> w0 = m_base[0]->GetW();
    Array<OneD, const NekDouble> w1 = m_base[1]->GetW();
 
    return StdExpansion2D::Integral(inarray, w0, w1);
}
 
/////////////////////////////
// Differentiation Methods //
/////////////////////////////
 
/** \brief Calculate the derivative of the physical points
 *
 *  For quadrilateral region can use the Tensor_Deriv function
 *  defined under StdExpansion.
 */
 
void StdQuadExp::v_PhysDeriv(const Array<OneD, const NekDouble> &inarray,
                             Array<OneD, NekDouble> &out_d0,
                             Array<OneD, NekDouble> &out_d1,
                             [[maybe_unused]] Array<OneD, NekDouble> &out_d2)
{
    PhysTensorDeriv(inarray, out_d0, out_d1);
}
 
void StdQuadExp::v_PhysDeriv(const int dir,
                             const Array<OneD, const NekDouble> &inarray,
                             Array<OneD, NekDouble> &outarray)
{
    switch (dir)
    {
        case 0:
        {
            PhysTensorDeriv(inarray, outarray, NullNekDouble1DArray);
        }
        break;
        case 1:
        {
            PhysTensorDeriv(inarray, NullNekDouble1DArray, outarray);
        }
        break;
        default:
        {
            ASSERTL1(false, "input dir is out of range");
        }
        break;
    }
}
 
void StdQuadExp::v_StdPhysDeriv(const Array<OneD, const NekDouble> &inarray,
                                Array<OneD, NekDouble> &out_d0,
                                Array<OneD, NekDouble> &out_d1,
                                [[maybe_unused]] Array<OneD, NekDouble> &out_d2)
{
    StdQuadExp::v_PhysDeriv(inarray, out_d0, out_d1);
}
 
////////////////
// Transforms //
////////////////
 
void StdQuadExp::v_BwdTrans(const Array<OneD, const NekDouble> &inarray,
                            Array<OneD, NekDouble> &outarray)
{
    if (m_base[0]->Collocation() && m_base[1]->Collocation())
    {
        Vmath::Vcopy(m_base[0]->GetNumPoints() * m_base[1]->GetNumPoints(),
                     inarray, 1, outarray, 1);
    }
    else
    {
        StdQuadExp::v_BwdTrans_SumFac(inarray, outarray);
    }
}
 
void StdQuadExp::v_BwdTrans_SumFac(const Array<OneD, const NekDouble> &inarray,
                                   Array<OneD, NekDouble> &outarray)
{
    Array<OneD, NekDouble> wsp(m_base[0]->GetNumPoints() *
                               m_base[1]->GetNumModes());
 
    BwdTrans_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetBdata(), inarray,
                          outarray, wsp, true, true);
}
 
// The arguments doCheckCollDir0 and doCheckCollDir1 allow you to specify
// whether to check if the basis has collocation properties (i.e. for the
// classical spectral element basis, In this case the 1D 'B' matrix is equal to
// the identity matrix which can be exploited to speed up the calculations).
// However, as this routine also allows to pass the matrix 'DB' (derivative of
// the basis), the collocation property cannot always be used. Therefor follow
// this rule: if base0 == m_base[0]->GetBdata() --> set doCheckCollDir0 == true;
//    base1 == m_base[1]->GetBdata() --> set doCheckCollDir1 == true;
//    base0 == m_base[0]->GetDbdata() --> set doCheckCollDir0 == false;
//    base1 == m_base[1]->GetDbdata() --> set doCheckCollDir1 == false;
void StdQuadExp::v_BwdTrans_SumFacKernel(
    const Array<OneD, const NekDouble> &base0,
    const Array<OneD, const NekDouble> &base1,
    const Array<OneD, const NekDouble> &inarray,
    Array<OneD, NekDouble> &outarray, Array<OneD, NekDouble> &wsp,
    bool doCheckCollDir0, bool doCheckCollDir1)
{
    int nquad0  = m_base[0]->GetNumPoints();
    int nquad1  = m_base[1]->GetNumPoints();
    int nmodes0 = m_base[0]->GetNumModes();
    int nmodes1 = m_base[1]->GetNumModes();
 
    bool colldir0 = doCheckCollDir0 ? (m_base[0]->Collocation()) : false;
    bool colldir1 = doCheckCollDir1 ? (m_base[1]->Collocation()) : false;
 
    if (colldir0 && colldir1)
    {
        Vmath::Vcopy(m_ncoeffs, inarray.data(), 1, outarray.data(), 1);
    }
    else if (colldir0)
    {
        Blas::Dgemm('N', 'T', nquad0, nquad1, nmodes1, 1.0, &inarray[0], nquad0,
                    base1.data(), nquad1, 0.0, &outarray[0], nquad0);
    }
    else if (colldir1)
    {
        Blas::Dgemm('N', 'N', nquad0, nmodes1, nmodes0, 1.0, base0.data(),
                    nquad0, &inarray[0], nmodes0, 0.0, &outarray[0], nquad0);
    }
    else
    {
        ASSERTL1(wsp.size() >= nquad0 * nmodes1,
                 "Workspace size is not sufficient");
 
        // Those two calls correpsond to the operation
        // out = B0*in*Transpose(B1);
        Blas::Dgemm('N', 'N', nquad0, nmodes1, nmodes0, 1.0, base0.data(),
                    nquad0, &inarray[0], nmodes0, 0.0, &wsp[0], nquad0);
        Blas::Dgemm('N', 'T', nquad0, nquad1, nmodes1, 1.0, &wsp[0], nquad0,
                    base1.data(), nquad1, 0.0, &outarray[0], nquad0);
    }
}
 
void StdQuadExp::v_FwdTrans(const Array<OneD, const NekDouble> &inarray,
                            Array<OneD, NekDouble> &outarray)
{
    if ((m_base[0]->Collocation()) && (m_base[1]->Collocation()))
    {
        Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
    }
    else
    {
        StdQuadExp::v_IProductWRTBase(inarray, outarray);
 
        // get Mass matrix inverse
        StdMatrixKey masskey(eInvMass, DetShapeType(), *this);
        DNekMatSharedPtr matsys = GetStdMatrix(masskey);
 
        // copy inarray in case inarray == outarray
        NekVector<NekDouble> in(m_ncoeffs, outarray, eCopy);
        NekVector<NekDouble> out(m_ncoeffs, outarray, eWrapper);
 
        out = (*matsys) * in;
    }
}
 
void StdQuadExp::v_FwdTransBndConstrained(
    const Array<OneD, const NekDouble> &inarray,
    Array<OneD, NekDouble> &outarray)
{
    if ((m_base[0]->Collocation()) && (m_base[1]->Collocation()))
    {
        Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
    }
    else
    {
        int i, j;
        int npoints[2] = {m_base[0]->GetNumPoints(), m_base[1]->GetNumPoints()};
        int nmodes[2]  = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes()};
 
        fill(outarray.data(), outarray.data() + m_ncoeffs, 0.0);
 
        Array<OneD, NekDouble> physEdge[4];
        Array<OneD, NekDouble> coeffEdge[4];
        for (i = 0; i < 4; i++)
        {
            physEdge[i]  = Array<OneD, NekDouble>(npoints[i % 2]);
            coeffEdge[i] = Array<OneD, NekDouble>(nmodes[i % 2]);
        }
 
        for (i = 0; i < npoints[0]; i++)
        {
            physEdge[0][i] = inarray[i];
            physEdge[2][i] = inarray[npoints[0] * npoints[1] - 1 - i];
        }
 
        for (i = 0; i < npoints[1]; i++)
        {
            physEdge[1][i] = inarray[npoints[0] - 1 + i * npoints[0]];
            physEdge[3][i] =
                inarray[(npoints[1] - 1) * npoints[0] - i * npoints[0]];
        }
 
        StdSegExpSharedPtr segexp[2] = {
            MemoryManager<StdRegions::StdSegExp>::AllocateSharedPtr(
                m_base[0]->GetBasisKey()),
            MemoryManager<StdRegions::StdSegExp>::AllocateSharedPtr(
                m_base[1]->GetBasisKey())};
 
        Array<OneD, unsigned int> mapArray;
        Array<OneD, int> signArray;
        NekDouble sign;
 
        for (i = 0; i < 4; i++)
        {
            segexp[i % 2]->FwdTransBndConstrained(physEdge[i], coeffEdge[i]);
 
            GetTraceToElementMap(i, mapArray, signArray);
            for (j = 0; j < nmodes[i % 2]; j++)
            {
                sign                  = (NekDouble)signArray[j];
                outarray[mapArray[j]] = sign * coeffEdge[i][j];
            }
        }
 
        Array<OneD, NekDouble> tmp0(m_ncoeffs);
        Array<OneD, NekDouble> tmp1(m_ncoeffs);
 
        StdMatrixKey masskey(eMass, DetShapeType(), *this);
        MassMatrixOp(outarray, tmp0, masskey);
        IProductWRTBase(inarray, tmp1);
 
        Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
 
        // get Mass matrix inverse (only of interior DOF)
        // use block (1,1) of the static condensed system
        // note: this block alreay contains the inverse matrix
        DNekMatSharedPtr matsys =
            (m_stdStaticCondMatrixManager[masskey])->GetBlock(1, 1);
 
        int nBoundaryDofs = NumBndryCoeffs();
        int nInteriorDofs = m_ncoeffs - nBoundaryDofs;
 
        Array<OneD, NekDouble> rhs(nInteriorDofs);
        Array<OneD, NekDouble> result(nInteriorDofs);
 
        GetInteriorMap(mapArray);
 
        for (i = 0; i < nInteriorDofs; i++)
        {
            rhs[i] = tmp1[mapArray[i]];
        }
 
        Blas::Dgemv('N', nInteriorDofs, nInteriorDofs, 1.0,
                    &(matsys->GetPtr())[0], nInteriorDofs, rhs.data(), 1, 0.0,
                    result.data(), 1);
 
        for (i = 0; i < nInteriorDofs; i++)
        {
            outarray[mapArray[i]] = result[i];
        }
    }
}
 
/////////////////////////////
// Inner Product Functions //
/////////////////////////////
 
/** \brief Calculate the inner product of inarray with respect to
 *  the basis B=base0*base1 and put into outarray
 *
 *  \f$
 *  \begin{array}{rcl}
 *  I_{pq} = (\phi_p \phi_q, u) & = & \sum_{i=0}^{nq_0}
 *  \sum_{j=0}^{nq_1}
 *  \phi_p(\xi_{0,i}) \phi_q(\xi_{1,j}) w^0_i w^1_j u(\xi_{0,i}
 *  \xi_{1,j}) \\
 *  & = & \sum_{i=0}^{nq_0} \phi_p(\xi_{0,i})
 *  \sum_{j=0}^{nq_1} \phi_q(\xi_{1,j}) \tilde{u}_{i,j}
 *  \end{array}
 *  \f$
 *
 *  where
 *
 *  \f$  \tilde{u}_{i,j} = w^0_i w^1_j u(\xi_{0,i},\xi_{1,j}) \f$
 *
 *  which can be implemented as
 *
 *  \f$  f_{qi} = \sum_{j=0}^{nq_1} \phi_q(\xi_{1,j})
 *  \tilde{u}_{i,j} = {\bf B_1 U}  \f$
 *  \f$  I_{pq} = \sum_{i=0}^{nq_0} \phi_p(\xi_{0,i}) f_{qi} =
 *  {\bf B_0 F}  \f$
 */
void StdQuadExp::v_IProductWRTBase(const Array<OneD, const NekDouble> &inarray,
                                   Array<OneD, NekDouble> &outarray)
{
    if (m_base[0]->Collocation() && m_base[1]->Collocation())
    {
        MultiplyByQuadratureMetric(inarray, outarray);
    }
    else
    {
        StdQuadExp::v_IProductWRTBase_SumFac(inarray, outarray);
    }
}
 
void StdQuadExp::v_IProductWRTBase_SumFac(
    const Array<OneD, const NekDouble> &inarray,
    Array<OneD, NekDouble> &outarray, bool multiplybyweights)
{
    int nquad0 = m_base[0]->GetNumPoints();
    int nquad1 = m_base[1]->GetNumPoints();
    int order0 = m_base[0]->GetNumModes();
 
    if (multiplybyweights)
    {
        Array<OneD, NekDouble> tmp(nquad0 * nquad1 + nquad1 * order0);
        Array<OneD, NekDouble> wsp(tmp + nquad0 * nquad1);
 
        // multiply by integration constants
        MultiplyByQuadratureMetric(inarray, tmp);
        IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
                                     m_base[1]->GetBdata(), tmp, outarray, wsp,
                                     true, true);
    }
    else
    {
        Array<OneD, NekDouble> wsp(nquad1 * order0);
        IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
                                     m_base[1]->GetBdata(), inarray, outarray,
                                     wsp, true, true);
    }
}
 
void StdQuadExp::v_IProductWRTDerivBase(
    const int dir, const Array<OneD, const NekDouble> &inarray,
    Array<OneD, NekDouble> &outarray)
{
    v_IProductWRTDerivBase_SumFac(dir, inarray, outarray);
}
 
void StdQuadExp::v_IProductWRTDerivBase_SumFac(
    const int dir, const Array<OneD, const NekDouble> &inarray,
    Array<OneD, NekDouble> &outarray)
{
    ASSERTL0((dir == 0) || (dir == 1), "input dir is out of range");
 
    int nquad0 = m_base[0]->GetNumPoints();
    int nquad1 = m_base[1]->GetNumPoints();
    int nqtot  = nquad0 * nquad1;
    int order0 = m_base[0]->GetNumModes();
 
    Array<OneD, NekDouble> tmp(nqtot + nquad1 * order0);
    Array<OneD, NekDouble> wsp(tmp + nqtot);
 
    // multiply by integration constants
    MultiplyByQuadratureMetric(inarray, tmp);
 
    if (dir) // dir == 1
    {
        IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
                                     m_base[1]->GetDbdata(), tmp, outarray, wsp,
                                     true, false);
    }
    else // dir == 0
    {
        IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
                                     m_base[1]->GetBdata(), tmp, outarray, wsp,
                                     false, true);
    }
}
 
// the arguments doCheckCollDir0 and doCheckCollDir1 allow you to specify
// whether to check if the basis has collocation properties (i.e. for the
// classical spectral element basis, In this case the 1D 'B' matrix is equal to
// the identity matrix which can be exploited to speed up the calculations).
// However, as this routine also allows to pass the matrix 'DB' (derivative of
// the basis), the collocation property cannot always be used. Therefor follow
// this rule: if base0 == m_base[0]->GetBdata() --> set doCheckCollDir0 == true;
//    base1 == m_base[1]->GetBdata() --> set doCheckCollDir1 == true;
//    base0 == m_base[0]->GetDbdata() --> set doCheckCollDir0 == false;
//    base1 == m_base[1]->GetDbdata() --> set doCheckCollDir1 == false;
void StdQuadExp::v_IProductWRTBase_SumFacKernel(
    const Array<OneD, const NekDouble> &base0,
    const Array<OneD, const NekDouble> &base1,
    const Array<OneD, const NekDouble> &inarray,
    Array<OneD, NekDouble> &outarray, Array<OneD, NekDouble> &wsp,
    bool doCheckCollDir0, bool doCheckCollDir1)
{
    int nquad0  = m_base[0]->GetNumPoints();
    int nquad1  = m_base[1]->GetNumPoints();
    int nmodes0 = m_base[0]->GetNumModes();
    int nmodes1 = m_base[1]->GetNumModes();
 
    bool colldir0 = doCheckCollDir0 ? (m_base[0]->Collocation()) : false;
    bool colldir1 = doCheckCollDir1 ? (m_base[1]->Collocation()) : false;
 
    if (colldir0 && colldir1)
    {
        Vmath::Vcopy(m_ncoeffs, inarray.data(), 1, outarray.data(), 1);
    }
    else if (colldir0)
    {
        Blas::Dgemm('N', 'N', nmodes0, nmodes1, nquad1, 1.0, inarray.data(),
                    nmodes0, base1.data(), nquad1, 0.0, outarray.data(),
                    nmodes0);
    }
    else if (colldir1)
    {
        Blas::Dgemm('T', 'N', nmodes0, nquad1, nquad0, 1.0, base0.data(),
                    nquad0, inarray.data(), nquad0, 0.0, outarray.data(),
                    nmodes0);
    }
    else
    {
        ASSERTL1(wsp.size() >= nquad1 * nmodes0,
                 "Workspace size is not sufficient");
 
#if 1
        Blas::Dgemm('T', 'N', nmodes0, nquad1, nquad0, 1.0, base0.data(),
                    nquad0, inarray.data(), nquad0, 0.0, wsp.data(), nmodes0);
 
#else
        for (int i = 0; i < nmodes0; ++i)
        {
            for (int j = 0; j < nquad1; ++j)
            {
                wsp[j * nmodes0 + i] =
                    Blas::Ddot(nquad0, base0.data() + i * nquad0, 1,
                               inarray.data() + j * nquad0, 1);
            }
        }
#endif
        Blas::Dgemm('N', 'N', nmodes0, nmodes1, nquad1, 1.0, wsp.data(),
                    nmodes0, base1.data(), nquad1, 0.0, outarray.data(),
                    nmodes0);
    }
}
 
//////////////////////////
// Evaluation functions //
//////////////////////////
 
void StdQuadExp::v_LocCoordToLocCollapsed(
    const Array<OneD, const NekDouble> &xi, Array<OneD, NekDouble> &eta)
{
    eta[0] = xi[0];
    eta[1] = xi[1];
}
 
void StdQuadExp::v_LocCollapsedToLocCoord(
    const Array<OneD, const NekDouble> &eta, Array<OneD, NekDouble> &xi)
{
    xi[0] = eta[0];
    xi[1] = eta[1];
}
 
/** \brief Fill outarray with mode \a mode of expansion
 *
 *  Note for quadrilateral expansions _base[0] (i.e. p)  modes run
 *  fastest
 */
 
void StdQuadExp::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
{
    int i;
    int nquad0                         = m_base[0]->GetNumPoints();
    int nquad1                         = m_base[1]->GetNumPoints();
    Array<OneD, const NekDouble> base0 = m_base[0]->GetBdata();
    Array<OneD, const NekDouble> base1 = m_base[1]->GetBdata();
    int btmp0                          = m_base[0]->GetNumModes();
    int mode0                          = mode % btmp0;
    int mode1                          = mode / btmp0;
 
    ASSERTL2(mode1 == (int)floor((1.0 * mode) / btmp0),
             "Integer Truncation not Equiv to Floor");
 
    ASSERTL2(m_ncoeffs > mode,
             "calling argument mode is larger than total expansion order");
 
    for (i = 0; i < nquad1; ++i)
    {
        Vmath::Vcopy(nquad0, (NekDouble *)(base0.data() + mode0 * nquad0), 1,
                     &outarray[0] + i * nquad0, 1);
    }
 
    for (i = 0; i < nquad0; ++i)
    {
        Vmath::Vmul(nquad1, (NekDouble *)(base1.data() + mode1 * nquad1), 1,
                    &outarray[0] + i, nquad0, &outarray[0] + i, nquad0);
    }
}
 
//////////////////////
// Helper functions //
//////////////////////
 
int StdQuadExp::v_GetNverts() const
{
    return 4;
}
 
int StdQuadExp::v_GetNtraces() const
{
    return 4;
}
 
int StdQuadExp::v_GetTraceNcoeffs(const int i) const
{
    ASSERTL2((i >= 0) && (i <= 3), "edge id is out of range");
 
    if ((i == 0) || (i == 2))
    {
        return GetBasisNumModes(0);
    }
    else
    {
        return GetBasisNumModes(1);
    }
}
 
int StdQuadExp::v_GetTraceIntNcoeffs(const int i) const
{
    ASSERTL2((i >= 0) && (i <= 4), "edge id is out of range");
    if ((i == 0) || (i == 2))
    {
        return GetBasisNumModes(0) - 2;
    }
    else
    {
        return GetBasisNumModes(1) - 2;
    }
}
 
int StdQuadExp::v_GetTraceNumPoints(const int i) const
{
    ASSERTL2((i >= 0) && (i <= 3), "edge id is out of range");
 
    if ((i == 0) || (i == 2))
    {
        return GetNumPoints(0);
    }
    else
    {
        return GetNumPoints(1);
    }
}
 
const LibUtilities::BasisKey StdQuadExp::v_GetTraceBasisKey(
    const int i, [[maybe_unused]] const int j,
    [[maybe_unused]] bool UseGLL) const
{
    ASSERTL2((i >= 0) && (i <= 3), "edge id is out of range");
 
    if ((i == 0) || (i == 2))
    {
        switch (GetBasis(0)->GetBasisType())
        {
            case LibUtilities::eOrtho_A:
            {
                return LibUtilities::BasisKey(LibUtilities::eGLL_Lagrange,
                                              GetBasis(0)->GetNumModes(),
                                              GetBasis(0)->GetPointsKey());
            }
            break;
            default:
            {
                return GetBasis(0)->GetBasisKey();
            }
        }
    }
    else
    {
        switch (GetBasis(1)->GetBasisType())
        {
            case LibUtilities::eOrtho_A:
            {
                return LibUtilities::BasisKey(LibUtilities::eGLL_Lagrange,
                                              GetBasis(1)->GetNumModes(),
                                              GetBasis(1)->GetPointsKey());
            }
            break;
            default:
            {
                return GetBasis(1)->GetBasisKey();
            }
        }
    }
}
 
LibUtilities::ShapeType StdQuadExp::v_DetShapeType() const
{
    return LibUtilities::eQuadrilateral;
}
 
int StdQuadExp::v_NumBndryCoeffs() const
{
    ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
                 GetBasisType(0) == LibUtilities::eGLL_Lagrange ||
                 GetBasisType(0) == LibUtilities::eGauss_Lagrange,
             "BasisType is not a boundary interior form");
    ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
                 GetBasisType(1) == LibUtilities::eGLL_Lagrange ||
                 GetBasisType(0) == LibUtilities::eGauss_Lagrange,
             "BasisType is not a boundary interior form");
 
    return 4 + 2 * (GetBasisNumModes(0) - 2) + 2 * (GetBasisNumModes(1) - 2);
}
 
int StdQuadExp::v_NumDGBndryCoeffs() const
{
    ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
                 GetBasisType(0) == LibUtilities::eGLL_Lagrange ||
                 GetBasisType(0) == LibUtilities::eGauss_Lagrange,
             "BasisType is not a boundary interior form");
    ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
                 GetBasisType(1) == LibUtilities::eGLL_Lagrange ||
                 GetBasisType(0) == LibUtilities::eGauss_Lagrange,
             "BasisType is not a boundary interior form");
 
    return 2 * GetBasisNumModes(0) + 2 * GetBasisNumModes(1);
}
 
int StdQuadExp::v_CalcNumberOfCoefficients(
    const std::vector<unsigned int> &nummodes, int &modes_offset)
{
    int nmodes = nummodes[modes_offset] * nummodes[modes_offset + 1];
    modes_offset += 2;
 
    return nmodes;
}
 
bool StdQuadExp::v_IsBoundaryInteriorExpansion() const
{
    bool returnval = false;
 
    if ((m_base[0]->GetBasisType() == LibUtilities::eModified_A) ||
        (m_base[0]->GetBasisType() == LibUtilities::eGLL_Lagrange))
    {
        if ((m_base[1]->GetBasisType() == LibUtilities::eModified_A) ||
            (m_base[1]->GetBasisType() == LibUtilities::eGLL_Lagrange))
        {
            returnval = true;
        }
    }
 
    return returnval;
}
 
void StdQuadExp::v_GetCoords(Array<OneD, NekDouble> &coords_0,
                             Array<OneD, NekDouble> &coords_1,
                             [[maybe_unused]] Array<OneD, NekDouble> &coords_2)
{
    Array<OneD, const NekDouble> z0 = m_base[0]->GetZ();
    Array<OneD, const NekDouble> z1 = m_base[1]->GetZ();
    int nq0                         = GetNumPoints(0);
    int nq1                         = GetNumPoints(1);
    int i;
 
    for (i = 0; i < nq1; ++i)
    {
        Blas::Dcopy(nq0, z0.data(), 1, &coords_0[0] + i * nq0, 1);
        Vmath::Fill(nq0, z1[i], &coords_1[0] + i * nq0, 1);
    }
}
 
/**
 * @brief This function evaluates the basis function mode @p mode at a
 * point @p coords of the domain.
 *
 * This function uses barycentric interpolation with the tensor
 * product separation of the basis function to improve performance.
 *
 * @param coord   The coordinate inside the standard region.
 * @param mode    The mode number to be evaluated.
 *
 * @return The value of the basis function @p mode at @p coords.
 */
NekDouble StdQuadExp::v_PhysEvaluateBasis(
    const Array<OneD, const NekDouble> &coords, int mode)
{
    ASSERTL2(coords[0] > -1 - NekConstants::kNekZeroTol, "coord[0] < -1");
    ASSERTL2(coords[0] < 1 + NekConstants::kNekZeroTol, "coord[0] >  1");
    ASSERTL2(coords[1] > -1 - NekConstants::kNekZeroTol, "coord[1] < -1");
    ASSERTL2(coords[1] < 1 + NekConstants::kNekZeroTol, "coord[1] >  1");
 
    const int nm0 = m_base[0]->GetNumModes();
    const int nm1 = m_base[1]->GetNumModes();
 
    return StdExpansion::BaryEvaluateBasis<0>(coords[0], mode % nm1) *
           StdExpansion::BaryEvaluateBasis<1>(coords[1], mode / nm0);
}
 
NekDouble StdQuadExp::v_PhysEvalFirstDeriv(
    const Array<OneD, NekDouble> &coord,
    const Array<OneD, const NekDouble> &inarray,
    std::array<NekDouble, 3> &firstOrderDerivs)
{
    return BaryTensorDeriv(coord, inarray, firstOrderDerivs);
}
 
//////////////
// Mappings //
//////////////
 
void StdQuadExp::v_GetBoundaryMap(Array<OneD, unsigned int> &outarray)
{
    int i;
    int cnt = 0;
    int nummodes0, nummodes1;
    int value1 = 0, value2 = 0;
    if (outarray.size() != NumBndryCoeffs())
    {
        outarray = Array<OneD, unsigned int>(NumBndryCoeffs());
    }
 
    nummodes0 = m_base[0]->GetNumModes();
    nummodes1 = m_base[1]->GetNumModes();
 
    const LibUtilities::BasisType Btype0 = GetBasisType(0);
    const LibUtilities::BasisType Btype1 = GetBasisType(1);
 
    switch (Btype1)
    {
        case LibUtilities::eGLL_Lagrange:
        case LibUtilities::eGauss_Lagrange:
            value1 = nummodes0;
            break;
        case LibUtilities::eModified_A:
            value1 = 2 * nummodes0;
            break;
        default:
            ASSERTL0(0, "Mapping array is not defined for this expansion");
            break;
    }
 
    for (i = 0; i < value1; i++)
    {
        outarray[i] = i;
    }
    cnt = value1;
 
    switch (Btype0)
    {
        case LibUtilities::eGLL_Lagrange:
        case LibUtilities::eGauss_Lagrange:
            value2 = value1 + nummodes0 - 1;
            break;
        case LibUtilities::eModified_A:
            value2 = value1 + 1;
            break;
        default:
            ASSERTL0(0, "Mapping array is not defined for this expansion");
            break;
    }
 
    for (i = 0; i < nummodes1 - 2; i++)
    {
        outarray[cnt++] = value1 + i * nummodes0;
        outarray[cnt++] = value2 + i * nummodes0;
    }
 
    if (Btype1 == LibUtilities::eGLL_Lagrange ||
        Btype1 == LibUtilities::eGauss_Lagrange)
    {
        for (i = nummodes0 * (nummodes1 - 1); i < GetNcoeffs(); i++)
        {
            outarray[cnt++] = i;
        }
    }
}
 
void StdQuadExp::v_GetInteriorMap(Array<OneD, unsigned int> &outarray)
{
    int i, j;
    int cnt = 0;
    int nummodes0, nummodes1;
    int startvalue = 0;
    if (outarray.size() != GetNcoeffs() - NumBndryCoeffs())
    {
        outarray = Array<OneD, unsigned int>(GetNcoeffs() - NumBndryCoeffs());
    }
 
    nummodes0 = m_base[0]->GetNumModes();
    nummodes1 = m_base[1]->GetNumModes();
 
    const LibUtilities::BasisType Btype0 = GetBasisType(0);
    const LibUtilities::BasisType Btype1 = GetBasisType(1);
 
    switch (Btype1)
    {
        case LibUtilities::eGLL_Lagrange:
            startvalue = nummodes0;
            break;
        case LibUtilities::eModified_A:
            startvalue = 2 * nummodes0;
            break;
        default:
            ASSERTL0(0, "Mapping array is not defined for this expansion");
            break;
    }
 
    switch (Btype0)
    {
        case LibUtilities::eGLL_Lagrange:
            startvalue++;
            break;
        case LibUtilities::eModified_A:
            startvalue += 2;
            break;
        default:
            ASSERTL0(0, "Mapping array is not defined for this expansion");
            break;
    }
 
    for (i = 0; i < nummodes1 - 2; i++)
    {
        for (j = 0; j < nummodes0 - 2; j++)
        {
            outarray[cnt++] = startvalue + j;
        }
        startvalue += nummodes0;
    }
}
 
int StdQuadExp::v_GetVertexMap(int localVertexId, bool useCoeffPacking)
{
    int localDOF = 0;
 
    if (useCoeffPacking == true)
    {
        switch (localVertexId)
        {
            case 0:
            {
                localDOF = 0;
            }
            break;
            case 1:
            {
                if (m_base[0]->GetBasisType() == LibUtilities::eGLL_Lagrange)
                {
                    localDOF = m_base[0]->GetNumModes() - 1;
                }
                else
                {
                    localDOF = 1;
                }
            }
            break;
            case 2:
            {
                if (m_base[0]->GetBasisType() == LibUtilities::eGLL_Lagrange)
                {
                    localDOF = m_base[0]->GetNumModes() *
                               (m_base[1]->GetNumModes() - 1);
                }
                else
                {
                    localDOF = m_base[0]->GetNumModes();
                }
            }
            break;
            case 3:
            {
                if (m_base[0]->GetBasisType() == LibUtilities::eGLL_Lagrange)
                {
                    localDOF =
                        m_base[0]->GetNumModes() * m_base[1]->GetNumModes() - 1;
                }
                else
                {
                    localDOF = m_base[0]->GetNumModes() + 1;
                }
            }
            break;
            default:
                ASSERTL0(false, "eid must be between 0 and 3");
                break;
        }
    }
    else
    {
        switch (localVertexId)
        {
            case 0:
            {
                localDOF = 0;
            }
            break;
            case 1:
            {
                if (m_base[0]->GetBasisType() == LibUtilities::eGLL_Lagrange)
                {
                    localDOF = m_base[0]->GetNumModes() - 1;
                }
                else
                {
                    localDOF = 1;
                }
            }
            break;
            case 2:
            {
                if (m_base[0]->GetBasisType() == LibUtilities::eGLL_Lagrange)
                {
                    localDOF =
                        m_base[0]->GetNumModes() * m_base[1]->GetNumModes() - 1;
                }
                else
                {
                    localDOF = m_base[0]->GetNumModes() + 1;
                }
            }
            break;
            case 3:
            {
                if (m_base[0]->GetBasisType() == LibUtilities::eGLL_Lagrange)
                {
                    localDOF = m_base[0]->GetNumModes() *
                               (m_base[1]->GetNumModes() - 1);
                }
                else
                {
                    localDOF = m_base[0]->GetNumModes();
                }
            }
            break;
            default:
                ASSERTL0(false, "eid must be between 0 and 3");
                break;
        }
    }
    return localDOF;
}
 
/**  \brief Get the map of the coefficient location to teh
 *    local trace coefficients
 */
 
void StdQuadExp::v_GetTraceCoeffMap(const unsigned int traceid,
                                    Array<OneD, unsigned int> &maparray)
{
    ASSERTL1(traceid < 4, "traceid must be between 0 and 3");
 
    unsigned int i;
    unsigned int order0   = m_base[0]->GetNumModes();
    unsigned int order1   = m_base[1]->GetNumModes();
    unsigned int numModes = (traceid % 2) ? order1 : order0;
 
    if (maparray.size() != numModes)
    {
        maparray = Array<OneD, unsigned int>(numModes);
    }
 
    const LibUtilities::BasisType bType = GetBasisType(traceid % 2);
 
    if (bType == LibUtilities::eModified_A)
    {
        switch (traceid)
        {
            case 0:
            {
                for (i = 0; i < numModes; i++)
                {
                    maparray[i] = i;
                }
            }
            break;
            case 1:
            {
                for (i = 0; i < numModes; i++)
                {
                    maparray[i] = i * order0 + 1;
                }
            }
            break;
            case 2:
            {
                for (i = 0; i < numModes; i++)
                {
                    maparray[i] = order0 + i;
                }
            }
            break;
            case 3:
            {
                for (i = 0; i < numModes; i++)
                {
                    maparray[i] = i * order0;
                }
            }
            break;
            default:
                break;
        }
    }
    else if (bType == LibUtilities::eGLL_Lagrange ||
             bType == LibUtilities::eGauss_Lagrange)
    {
        switch (traceid)
        {
            case 0:
            {
                for (i = 0; i < numModes; i++)
                {
                    maparray[i] = i;
                }
            }
            break;
            case 1:
            {
                for (i = 0; i < numModes; i++)
                {
                    maparray[i] = (i + 1) * order0 - 1;
                }
            }
            break;
            case 2:
            {
                for (i = 0; i < numModes; i++)
                {
                    maparray[i] = order0 * (order1 - 1) + i;
                }
            }
            break;
            case 3:
            {
                for (i = 0; i < numModes; i++)
                {
                    maparray[i] = order0 * i;
                }
            }
            break;
            default:
                break;
        }
    }
    else
    {
        ASSERTL0(false, "Mapping not defined for this type of basis");
    }
}
 
void StdQuadExp::v_GetTraceInteriorToElementMap(
    const int eid, Array<OneD, unsigned int> &maparray,
    Array<OneD, int> &signarray, const Orientation edgeOrient)
{
    int i;
    const int nummodes0                 = m_base[0]->GetNumModes();
    const int nummodes1                 = m_base[1]->GetNumModes();
    const int nEdgeIntCoeffs            = GetTraceNcoeffs(eid) - 2;
    const LibUtilities::BasisType bType = GetBasisType(eid % 2);
 
    if (maparray.size() != nEdgeIntCoeffs)
    {
        maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
    }
 
    if (signarray.size() != nEdgeIntCoeffs)
    {
        signarray = Array<OneD, int>(nEdgeIntCoeffs, 1);
    }
    else
    {
        fill(signarray.data(), signarray.data() + nEdgeIntCoeffs, 1);
    }
 
    if (bType == LibUtilities::eModified_A)
    {
        switch (eid)
        {
            case 0:
            {
                for (i = 0; i < nEdgeIntCoeffs; i++)
                {
                    maparray[i] = i + 2;
                }
            }
            break;
            case 1:
            {
                for (i = 0; i < nEdgeIntCoeffs; i++)
                {
                    maparray[i] = (i + 2) * nummodes0 + 1;
                }
            }
            break;
            case 2:
            {
                for (i = 0; i < nEdgeIntCoeffs; i++)
                {
                    maparray[i] = nummodes0 + i + 2;
                }
            }
            break;
            case 3:
            {
                for (i = 0; i < nEdgeIntCoeffs; i++)
                {
                    maparray[i] = (i + 2) * nummodes0;
                }
            }
            break;
            default:
                ASSERTL0(false, "eid must be between 0 and 3");
                break;
        }
 
        if (edgeOrient == eBackwards)
        {
            for (i = 1; i < nEdgeIntCoeffs; i += 2)
            {
                signarray[i] = -1;
            }
        }
    }
    else if (bType == LibUtilities::eGLL_Lagrange)
    {
        switch (eid)
        {
            case 0:
            {
                for (i = 0; i < nEdgeIntCoeffs; i++)
                {
                    maparray[i] = i + 1;
                }
            }
            break;
            case 1:
            {
                for (i = 0; i < nEdgeIntCoeffs; i++)
                {
                    maparray[i] = (i + 2) * nummodes0 - 1;
                }
            }
            break;
            case 2:
            {
                for (i = 0; i < nEdgeIntCoeffs; i++)
                {
                    maparray[i] = nummodes0 * (nummodes1 - 1) + i + 1;
                }
            }
            break;
            case 3:
            {
                for (i = 0; i < nEdgeIntCoeffs; i++)
                {
                    maparray[i] = nummodes0 * (i + 1);
                }
            }
            break;
            default:
                ASSERTL0(false, "eid must be between 0 and 3");
                break;
        }
        if (edgeOrient == eBackwards)
        {
            reverse(maparray.data(), maparray.data() + nEdgeIntCoeffs);
        }
    }
    else
    {
        ASSERTL0(false, "Mapping not defined for this type of basis");
    }
}
 
///////////////////////
// Wrapper Functions //
///////////////////////
 
DNekMatSharedPtr StdQuadExp::v_GenMatrix(const StdMatrixKey &mkey)
{
    int i, j;
    int order0       = GetBasisNumModes(0);
    int order1       = GetBasisNumModes(1);
    MatrixType mtype = mkey.GetMatrixType();
 
    DNekMatSharedPtr Mat;
 
    switch (mtype)
    {
        case ePhysInterpToEquiSpaced:
        {
            int nq0 = m_base[0]->GetNumPoints();
            int nq1 = m_base[1]->GetNumPoints();
            int nq;
 
            // take definition from key
            if (mkey.ConstFactorExists(eFactorConst))
            {
                nq = (int)mkey.GetConstFactor(eFactorConst);
            }
            else
            {
                nq = max(nq0, nq1);
            }
 
            int neq =
                LibUtilities::StdQuadData::getNumberOfCoefficients(nq, nq);
            Array<OneD, Array<OneD, NekDouble>> coords(neq);
            Array<OneD, NekDouble> coll(2);
            Array<OneD, DNekMatSharedPtr> I(2);
            Array<OneD, NekDouble> tmp(nq0);
 
            Mat     = MemoryManager<DNekMat>::AllocateSharedPtr(neq, nq0 * nq1);
            int cnt = 0;
 
            for (i = 0; i < nq; ++i)
            {
                for (j = 0; j < nq; ++j, ++cnt)
                {
                    coords[cnt]    = Array<OneD, NekDouble>(2);
                    coords[cnt][0] = -1.0 + 2 * j / (NekDouble)(nq - 1);
                    coords[cnt][1] = -1.0 + 2 * i / (NekDouble)(nq - 1);
                }
            }
 
            for (i = 0; i < neq; ++i)
            {
                LocCoordToLocCollapsed(coords[i], coll);
 
                I[0] = m_base[0]->GetI(coll);
                I[1] = m_base[1]->GetI(coll + 1);
 
                // interpolate first coordinate direction
                for (j = 0; j < nq1; ++j)
                {
                    NekDouble fac = (I[1]->GetPtr())[j];
                    Vmath::Smul(nq0, fac, I[0]->GetPtr(), 1, tmp, 1);
 
                    Vmath::Vcopy(nq0, &tmp[0], 1,
                                 Mat->GetRawPtr() + j * nq0 * neq + i, neq);
                }
            }
            break;
        }
        case eMass:
        {
            Mat = StdExpansion::CreateGeneralMatrix(mkey);
            // For Fourier basis set the imaginary component of mean mode
            // to have a unit diagonal component in mass matrix
            if (m_base[0]->GetBasisType() == LibUtilities::eFourier)
            {
                for (i = 0; i < order1; ++i)
                {
                    (*Mat)(order0 *i + 1, i * order0 + 1) = 1.0;
                }
            }
 
            if (m_base[1]->GetBasisType() == LibUtilities::eFourier)
            {
                for (i = 0; i < order0; ++i)
                {
                    (*Mat)(order0 + i, order0 + i) = 1.0;
                }
            }
            break;
        }
        case eFwdTrans:
        {
            Mat =
                MemoryManager<DNekMat>::AllocateSharedPtr(m_ncoeffs, m_ncoeffs);
            StdMatrixKey iprodkey(eIProductWRTBase, DetShapeType(), *this);
            DNekMat &Iprod = *GetStdMatrix(iprodkey);
            StdMatrixKey imasskey(eInvMass, DetShapeType(), *this);
            DNekMat &Imass = *GetStdMatrix(imasskey);
 
            (*Mat) = Imass * Iprod;
            break;
        }
        case eGaussDG:
        {
            ConstFactorMap factors = mkey.GetConstFactors();
 
            int edge    = (int)factors[StdRegions::eFactorGaussEdge];
            int dir     = (edge + 1) % 2;
            int nCoeffs = m_base[dir]->GetNumModes();
 
            const LibUtilities::PointsKey BS_p(
                nCoeffs, LibUtilities::eGaussGaussLegendre);
            const LibUtilities::BasisKey BS_k(LibUtilities::eGauss_Lagrange,
                                              nCoeffs, BS_p);
 
            Array<OneD, NekDouble> coords(1, 0.0);
            coords[0] = (edge == 0 || edge == 3) ? -1.0 : 1.0;
 
            LibUtilities::BasisSharedPtr basis =
                LibUtilities::BasisManager()[BS_k];
            DNekMatSharedPtr m_Ix = basis->GetI(coords);
 
            Mat = MemoryManager<DNekMat>::AllocateSharedPtr(1.0, nCoeffs);
            Vmath::Vcopy(nCoeffs, m_Ix->GetPtr(), 1, Mat->GetPtr(), 1);
            break;
        }
        default:
        {
            Mat = StdExpansion::CreateGeneralMatrix(mkey);
            break;
        }
    }
 
    return Mat;
}
 
DNekMatSharedPtr StdQuadExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
{
    return GenMatrix(mkey);
}
 
///////////////////////////////////
// Operator evaluation functions //
///////////////////////////////////
 
void StdQuadExp::v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
                                      const StdMatrixKey &mkey)
{
    // Generate an orthonogal expansion
    int qa       = m_base[0]->GetNumPoints();
    int qb       = m_base[1]->GetNumPoints();
    int nmodes_a = m_base[0]->GetNumModes();
    int nmodes_b = m_base[1]->GetNumModes();
    int nmodes   = min(nmodes_a, nmodes_b);
    // Declare orthogonal basis.
    LibUtilities::PointsKey pa(qa, m_base[0]->GetPointsType());
    LibUtilities::PointsKey pb(qb, m_base[1]->GetPointsType());
 
    LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, nmodes_a, pa);
    LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, nmodes_b, pb);
    StdQuadExp OrthoExp(Ba, Bb);
 
    // SVV parameters loaded from the .xml case file
    Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
 
    // project onto modal  space.
    OrthoExp.FwdTrans(array, orthocoeffs);
 
    if (mkey.ConstFactorExists(
            eFactorSVVPowerKerDiffCoeff)) // Rodrigo's power kernel
    {
        NekDouble cutoff = mkey.GetConstFactor(eFactorSVVCutoffRatio);
        NekDouble SvvDiffCoeff =
            mkey.GetConstFactor(eFactorSVVPowerKerDiffCoeff) *
            mkey.GetConstFactor(eFactorSVVDiffCoeff);
 
        for (int j = 0; j < nmodes_a; ++j)
        {
            for (int k = 0; k < nmodes_b; ++k)
            {
                // linear space but makes high modes very negative
                NekDouble fac = std::max(
                    pow((1.0 * j) / (nmodes_a - 1), cutoff * nmodes_a),
                    pow((1.0 * k) / (nmodes_b - 1), cutoff * nmodes_b));
 
                orthocoeffs[j * nmodes_b + k] *= SvvDiffCoeff * fac;
            }
        }
    }
    else if (mkey.ConstFactorExists(
                 eFactorSVVDGKerDiffCoeff)) // Rodrigo/mansoor's DG kernel
    {
        NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDGKerDiffCoeff) *
                                 mkey.GetConstFactor(eFactorSVVDiffCoeff);
        int max_ab = max(nmodes_a - kSVVDGFiltermodesmin,
                         nmodes_b - kSVVDGFiltermodesmin);
        max_ab     = max(max_ab, 0);
        max_ab     = min(max_ab, kSVVDGFiltermodesmax - kSVVDGFiltermodesmin);
 
        for (int j = 0; j < nmodes_a; ++j)
        {
            for (int k = 0; k < nmodes_b; ++k)
            {
                int maxjk = max(j, k);
                maxjk     = min(maxjk, kSVVDGFiltermodesmax - 1);
 
                orthocoeffs[j * nmodes_b + k] *=
                    SvvDiffCoeff * kSVVDGFilter[max_ab][maxjk];
            }
        }
    }
    else
    {
        NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDiffCoeff);
        // Exponential Kernel implementation
        int cutoff = (int)(mkey.GetConstFactor(eFactorSVVCutoffRatio) *
                           min(nmodes_a, nmodes_b));
 
        //------"New" Version August 22nd '13--------------------
        for (int j = 0; j < nmodes_a; ++j)
        {
            for (int k = 0; k < nmodes_b; ++k)
            {
                if (j + k >= cutoff) // to filter out only the "high-modes"
                {
                    orthocoeffs[j * nmodes_b + k] *=
                        (SvvDiffCoeff *
                         exp(-(j + k - nmodes) * (j + k - nmodes) /
                             ((NekDouble)((j + k - cutoff + 1) *
                                          (j + k - cutoff + 1)))));
                }
                else
                {
                    orthocoeffs[j * nmodes_b + k] *= 0.0;
                }
            }
        }
    }
 
    // backward transform to physical space
    OrthoExp.BwdTrans(orthocoeffs, array);
}
 
void StdQuadExp::v_ExponentialFilter(Array<OneD, NekDouble> &array,
                                     const NekDouble alpha,
                                     const NekDouble exponent,
                                     const NekDouble cutoff)
{
    // Generate an orthogonal expansion
    int qa      = m_base[0]->GetNumPoints();
    int qb      = m_base[1]->GetNumPoints();
    int nmodesA = m_base[0]->GetNumModes();
    int nmodesB = m_base[1]->GetNumModes();
    int P       = nmodesA - 1;
    int Q       = nmodesB - 1;
 
    // Declare orthogonal basis.
    LibUtilities::PointsKey pa(qa, m_base[0]->GetPointsType());
    LibUtilities::PointsKey pb(qb, m_base[1]->GetPointsType());
 
    LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, nmodesA, pa);
    LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, nmodesB, pb);
    StdQuadExp OrthoExp(Ba, Bb);
 
    // Cutoff
    int Pcut = cutoff * P;
    int Qcut = cutoff * Q;
 
    // Project onto orthogonal space.
    Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
    OrthoExp.FwdTrans(array, orthocoeffs);
 
    //
    NekDouble fac, fac1, fac2;
    for (int i = 0; i < nmodesA; ++i)
    {
        for (int j = 0; j < nmodesB; ++j)
        {
            // to filter out only the "high-modes"
            if (i > Pcut || j > Qcut)
            {
                fac1 = (NekDouble)(i - Pcut) / ((NekDouble)(P - Pcut));
                fac2 = (NekDouble)(j - Qcut) / ((NekDouble)(Q - Qcut));
                fac  = max(fac1, fac2);
                fac  = pow(fac, exponent);
                orthocoeffs[i * nmodesB + j] *= exp(-alpha * fac);
            }
        }
    }
 
    // backward transform to physical space
    OrthoExp.BwdTrans(orthocoeffs, array);
}
 
void StdQuadExp::v_ReduceOrderCoeffs(
    int numMin, const Array<OneD, const NekDouble> &inarray,
    Array<OneD, NekDouble> &outarray)
{
    int n_coeffs = inarray.size();
 
    Array<OneD, NekDouble> coeff(n_coeffs);
    Array<OneD, NekDouble> coeff_tmp(n_coeffs, 0.0);
    Array<OneD, NekDouble> tmp;
    Array<OneD, NekDouble> tmp2;
 
    int nmodes0 = m_base[0]->GetNumModes();
    int nmodes1 = m_base[1]->GetNumModes();
    int numMax  = nmodes0;
 
    Vmath::Vcopy(n_coeffs, inarray, 1, coeff_tmp, 1);
 
    const LibUtilities::PointsKey Pkey0(nmodes0,
                                        LibUtilities::eGaussLobattoLegendre);
    const LibUtilities::PointsKey Pkey1(nmodes1,
                                        LibUtilities::eGaussLobattoLegendre);
 
    LibUtilities::BasisKey b0(m_base[0]->GetBasisType(), nmodes0, Pkey0);
    LibUtilities::BasisKey b1(m_base[1]->GetBasisType(), nmodes1, Pkey1);
 
    LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A, nmodes0, Pkey0);
    LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_A, nmodes1, Pkey1);
 
    LibUtilities::InterpCoeff2D(b0, b1, coeff_tmp, bortho0, bortho1, coeff);
 
    Vmath::Zero(n_coeffs, coeff_tmp, 1);
 
    int cnt = 0;
    for (int i = 0; i < numMin + 1; ++i)
    {
        Vmath::Vcopy(numMin, tmp = coeff + cnt, 1, tmp2 = coeff_tmp + cnt, 1);
 
        cnt = i * numMax;
    }
 
    LibUtilities::InterpCoeff2D(bortho0, bortho1, coeff_tmp, b0, b1, outarray);
}
 
void StdQuadExp::v_MassMatrixOp(const Array<OneD, const NekDouble> &inarray,
                                Array<OneD, NekDouble> &outarray,
                                const StdMatrixKey &mkey)
{
    StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
}
 
void StdQuadExp::v_LaplacianMatrixOp(
    const Array<OneD, const NekDouble> &inarray,
    Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
{
    StdQuadExp::v_LaplacianMatrixOp_MatFree(inarray, outarray, mkey);
}
 
void StdQuadExp::v_LaplacianMatrixOp(
    const int k1, const int k2, const Array<OneD, const NekDouble> &inarray,
    Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
{
    StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
}
 
void StdQuadExp::v_WeakDerivMatrixOp(
    const int i, const Array<OneD, const NekDouble> &inarray,
    Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
{
    StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
}
 
void StdQuadExp::v_HelmholtzMatrixOp(
    const Array<OneD, const NekDouble> &inarray,
    Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
{
    StdQuadExp::v_HelmholtzMatrixOp_MatFree(inarray, outarray, mkey);
}
 
// up to here
void StdQuadExp::v_MultiplyByStdQuadratureMetric(
    const Array<OneD, const NekDouble> &inarray,
    Array<OneD, NekDouble> &outarray)
{
    int i;
    int nquad0 = m_base[0]->GetNumPoints();
    int nquad1 = m_base[1]->GetNumPoints();
 
    const Array<OneD, const NekDouble> &w0 = m_base[0]->GetW();
    const Array<OneD, const NekDouble> &w1 = m_base[1]->GetW();
 
    // multiply by integration constants
    for (i = 0; i < nquad1; ++i)
    {
        Vmath::Vmul(nquad0, inarray.data() + i * nquad0, 1, w0.data(), 1,
                    outarray.data() + i * nquad0, 1);
    }
 
    for (i = 0; i < nquad0; ++i)
    {
        Vmath::Vmul(nquad1, outarray.data() + i, nquad0, w1.data(), 1,
                    outarray.data() + i, nquad0);
    }
}
 
void StdQuadExp::v_GetSimplexEquiSpacedConnectivity(
    Array<OneD, int> &conn, [[maybe_unused]] bool standard)
{
    int np1 = m_base[0]->GetNumPoints();
    int np2 = m_base[1]->GetNumPoints();
    int np  = max(np1, np2);
 
    conn = Array<OneD, int>(6 * (np - 1) * (np - 1));
 
    int row   = 0;
    int rowp1 = 0;
    int cnt   = 0;
    for (int i = 0; i < np - 1; ++i)
    {
        rowp1 += np;
        for (int j = 0; j < np - 1; ++j)
        {
            conn[cnt++] = row + j;
            conn[cnt++] = row + j + 1;
            conn[cnt++] = rowp1 + j;
 
            conn[cnt++] = rowp1 + j + 1;
            conn[cnt++] = rowp1 + j;
            conn[cnt++] = row + j + 1;
        }
        row += np;
    }
}
 
} // namespace Nektar::StdRegions