StdTriExp.cpp

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File: StdTriExp.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: Triangle routines built upon StdExpansion2D
//
///////////////////////////////////////////////////////////////////////////////
 
#include <boost/core/ignore_unused.hpp>
 
#include <LibUtilities/Foundations/InterpCoeff.h>
#include <StdRegions/StdNodalTriExp.h>
#include <StdRegions/StdSegExp.h> // for StdSegExp, etc
#include <StdRegions/StdTriExp.h>
 
using namespace std;
 
namespace Nektar
{
namespace StdRegions
{
StdTriExp::StdTriExp()
{
}
 
StdTriExp::StdTriExp(const LibUtilities::BasisKey &Ba,
                     const LibUtilities::BasisKey &Bb)
    : StdExpansion(LibUtilities::StdTriData::getNumberOfCoefficients(
                       Ba.GetNumModes(), Bb.GetNumModes()),
                   2, Ba, Bb),
      StdExpansion2D(LibUtilities::StdTriData::getNumberOfCoefficients(
                         Ba.GetNumModes(), Bb.GetNumModes()),
                     Ba, Bb)
{
    ASSERTL0(Ba.GetNumModes() <= Bb.GetNumModes(),
             "order in 'a' direction is higher than order "
             "in 'b' direction");
}
 
StdTriExp::StdTriExp(const StdTriExp &T) : StdExpansion(T), StdExpansion2D(T)
{
}
 
// Destructor
StdTriExp::~StdTriExp()
{
}
 
//-------------------------------
// Integration Methods
//-------------------------------
NekDouble StdTriExp::v_Integral(const Array<OneD, const NekDouble> &inarray)
{
    int i;
    int nquad1 = m_base[1]->GetNumPoints();
    Array<OneD, NekDouble> w1_tmp(nquad1);
 
    Array<OneD, const NekDouble> w0 = m_base[0]->GetW();
    Array<OneD, const NekDouble> w1 = m_base[1]->GetW();
    Array<OneD, const NekDouble> z1 = m_base[1]->GetZ();
 
    switch (m_base[1]->GetPointsType())
    {
        case LibUtilities::eGaussRadauMAlpha1Beta0: // (0,1) Jacobi Inner
                                                    // product
        {
            Vmath::Smul(nquad1, 0.5, w1, 1, w1_tmp, 1);
            break;
        }
        default:
        {
            // include jacobian factor on whatever coordinates are defined.
            for (i = 0; i < nquad1; ++i)
            {
                w1_tmp[i] = 0.5 * (1 - z1[i]) * w1[i];
            }
            break;
        }
    }
 
    return StdExpansion2D::Integral(inarray, w0, w1_tmp);
}
 
//-----------------------------
// Differentiation Methods
//-----------------------------
 
/**
 * \brief Calculate the derivative of the physical points.
 *
 * \f$ \frac{\partial u}{\partial  x_1} =  \left .
 * \frac{2.0}{1-\eta_2} \frac{\partial u}{\partial d\eta_1}
 * \right |_{\eta_2}\f$
 *
 * \f$ \frac{\partial u}{\partial  x_2} =  \left .
 * \frac{1+\eta_1}{1-\eta_2} \frac{\partial u}{\partial d\eta_1}
 * \right |_{\eta_2}  + \left . \frac{\partial u}{\partial d\eta_2}
 * \right |_{\eta_1}  \f$
 */
void StdTriExp::v_PhysDeriv(const Array<OneD, const NekDouble> &inarray,
                            Array<OneD, NekDouble> &out_d0,
                            Array<OneD, NekDouble> &out_d1,
                            Array<OneD, NekDouble> &out_d2)
{
    boost::ignore_unused(out_d2);
 
    int i;
    int nquad0 = m_base[0]->GetNumPoints();
    int nquad1 = m_base[1]->GetNumPoints();
    Array<OneD, NekDouble> wsp(std::max(nquad0, nquad1));
 
    const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
    const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
 
    // set up geometric factor: 2/(1-z1)
    Vmath::Sadd(nquad1, -1.0, z1, 1, wsp, 1);
    Vmath::Sdiv(nquad1, -2.0, wsp, 1, wsp, 1);
 
    if (out_d0.size() > 0)
    {
        PhysTensorDeriv(inarray, out_d0, out_d1);
 
        for (i = 0; i < nquad1; ++i)
        {
            Blas::Dscal(nquad0, wsp[i], &out_d0[0] + i * nquad0, 1);
        }
 
        // if no d1 required do not need to calculate both deriv
        if (out_d1.size() > 0)
        {
            // set up geometric factor: (1+z0)/2
            Vmath::Sadd(nquad0, 1.0, z0, 1, wsp, 1);
            Vmath::Smul(nquad0, 0.5, wsp, 1, wsp, 1);
 
            for (i = 0; i < nquad1; ++i)
            {
                Vmath::Vvtvp(nquad0, &wsp[0], 1, &out_d0[0] + i * nquad0, 1,
                             &out_d1[0] + i * nquad0, 1,
                             &out_d1[0] + i * nquad0, 1);
            }
        }
    }
    else if (out_d1.size() > 0)
    {
        Array<OneD, NekDouble> diff0(nquad0 * nquad1);
        PhysTensorDeriv(inarray, diff0, out_d1);
 
        for (i = 0; i < nquad1; ++i)
        {
            Blas::Dscal(nquad0, wsp[i], &diff0[0] + i * nquad0, 1);
        }
 
        Vmath::Sadd(nquad0, 1.0, z0, 1, wsp, 1);
        Vmath::Smul(nquad0, 0.5, wsp, 1, wsp, 1);
 
        for (i = 0; i < nquad1; ++i)
        {
            Vmath::Vvtvp(nquad0, &wsp[0], 1, &diff0[0] + i * nquad0, 1,
                         &out_d1[0] + i * nquad0, 1, &out_d1[0] + i * nquad0,
                         1);
        }
    }
}
 
void StdTriExp::v_PhysDeriv(const int dir,
                            const Array<OneD, const NekDouble> &inarray,
                            Array<OneD, NekDouble> &outarray)
{
    switch (dir)
    {
        case 0:
        {
            v_PhysDeriv(inarray, outarray, NullNekDouble1DArray);
            break;
        }
        case 1:
        {
            v_PhysDeriv(inarray, NullNekDouble1DArray, outarray);
            break;
        }
        default:
        {
            ASSERTL1(false, "input dir is out of range");
            break;
        }
    }
}
 
void StdTriExp::v_StdPhysDeriv(const Array<OneD, const NekDouble> &inarray,
                               Array<OneD, NekDouble> &out_d0,
                               Array<OneD, NekDouble> &out_d1,
                               Array<OneD, NekDouble> &out_d2)
{
    boost::ignore_unused(out_d2);
    StdTriExp::v_PhysDeriv(inarray, out_d0, out_d1);
}
 
void StdTriExp::v_StdPhysDeriv(const int dir,
                               const Array<OneD, const NekDouble> &inarray,
                               Array<OneD, NekDouble> &outarray)
{
    StdTriExp::v_PhysDeriv(dir, inarray, outarray);
}
 
//---------------------------------------
// Transforms
//---------------------------------------
 
/**
 * \brief Backward tranform for triangular elements
 *
 * @note 'q' (base[1]) runs fastest in this element.
 */
void StdTriExp::v_BwdTrans(const Array<OneD, const NekDouble> &inarray,
                           Array<OneD, NekDouble> &outarray)
{
    v_BwdTrans_SumFac(inarray, outarray);
}
 
void StdTriExp::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);
}
 
void StdTriExp::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)
{
    boost::ignore_unused(doCheckCollDir0, doCheckCollDir1);
 
    int i;
    int mode;
    int nquad0  = m_base[0]->GetNumPoints();
    int nquad1  = m_base[1]->GetNumPoints();
    int nmodes0 = m_base[0]->GetNumModes();
    int nmodes1 = m_base[1]->GetNumModes();
 
    ASSERTL1(wsp.size() >= nquad0 * nmodes1,
             "Workspace size is not sufficient");
    ASSERTL2((m_base[1]->GetBasisType() != LibUtilities::eOrtho_B) ||
                 (m_base[1]->GetBasisType() != LibUtilities::eModified_B),
             "Basis[1] is not of general tensor type");
 
    for (i = mode = 0; i < nmodes0; ++i)
    {
        Blas::Dgemv('N', nquad1, nmodes1 - i, 1.0, base1.get() + mode * nquad1,
                    nquad1, &inarray[0] + mode, 1, 0.0, &wsp[0] + i * nquad1,
                    1);
        mode += nmodes1 - i;
    }
 
    // fix for modified basis by splitting top vertex mode
    if (m_base[0]->GetBasisType() == LibUtilities::eModified_A)
    {
        Blas::Daxpy(nquad1, inarray[1], base1.get() + nquad1, 1,
                    &wsp[0] + nquad1, 1);
    }
 
    Blas::Dgemm('N', 'T', nquad0, nquad1, nmodes0, 1.0, base0.get(), nquad0,
                &wsp[0], nquad1, 0.0, &outarray[0], nquad0);
}
 
void StdTriExp::v_FwdTrans(const Array<OneD, const NekDouble> &inarray,
                           Array<OneD, NekDouble> &outarray)
{
    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 StdTriExp::v_FwdTransBndConstrained(
    const Array<OneD, const NekDouble> &inarray,
    Array<OneD, NekDouble> &outarray)
{
    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.get(), outarray.get() + m_ncoeffs, 0.0);
 
    Array<OneD, NekDouble> physEdge[3];
    Array<OneD, NekDouble> coeffEdge[3];
    for (i = 0; i < 3; i++)
    {
        physEdge[i]  = Array<OneD, NekDouble>(npoints[i != 0]);
        coeffEdge[i] = Array<OneD, NekDouble>(nmodes[i != 0]);
    }
 
    for (i = 0; i < npoints[0]; i++)
    {
        physEdge[0][i] = inarray[i];
    }
 
    for (i = 0; i < npoints[1]; i++)
    {
        physEdge[1][i] = inarray[npoints[0] - 1 + i * npoints[0]];
        physEdge[2][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 < 3; i++)
    {
        segexp[i != 0]->FwdTransBndConstrained(physEdge[i], coeffEdge[i]);
 
        GetTraceToElementMap(i, mapArray, signArray);
        for (j = 0; j < nmodes[i != 0]; 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);
    v_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 = v_NumBndryCoeffs();
    int nInteriorDofs = m_ncoeffs - nBoundaryDofs;
 
    Array<OneD, NekDouble> rhs(nInteriorDofs);
    Array<OneD, NekDouble> result(nInteriorDofs);
 
    v_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.get(), 1, 0.0, result.get(), 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[p]*base1[pq] and put into outarray.
 *
 * \f$
 * \begin{array}{rcl}
 * I_{pq} = (\phi^A_q \phi^B_{pq}, u) &=&
 * \sum_{i=0}^{nq_0}\sum_{j=0}^{nq_1}
 * \phi^A_p(\eta_{0,i})\phi^B_{pq}(\eta_{1,j}) w^0_i w^1_j
 * u(\xi_{0,i} \xi_{1,j}) \\
 * & = & \sum_{i=0}^{nq_0} \phi^A_p(\eta_{0,i})
 * \sum_{j=0}^{nq_1} \phi^B_{pq}(\eta_{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_{pj} = \sum_{i=0}^{nq_0} \phi^A_p(\eta_{0,i})
 * \tilde{u}_{i,j}
 * \rightarrow {\bf B_1 U}  \f$
 * \f$  I_{pq} = \sum_{j=0}^{nq_1} \phi^B_{pq}(\eta_{1,j}) f_{pj}
 * \rightarrow {\bf B_2[p*skip] f[skip]}  \f$
 *
 * \b Recall: \f$ \eta_{1} = \frac{2(1+\xi_1)}{(1-\xi_2)}-1, \,
 * \eta_2 = \xi_2\f$
 *
 * \b Note: For the orthgonality of this expansion to be realised the
 * 'q' ordering must run fastest in contrast to the Quad and Hex
 * ordering where 'p' index runs fastest to be consistent with the
 * quadrature ordering.
 *
 * In the triangular space the i (i.e. \f$\eta_1\f$ direction)
 * ordering still runs fastest by convention.
 */
void StdTriExp::v_IProductWRTBase(const Array<OneD, const NekDouble> &inarray,
                                  Array<OneD, NekDouble> &outarray)
{
    StdTriExp::v_IProductWRTBase_SumFac(inarray, outarray);
}
 
void StdTriExp::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);
    }
    else
    {
        Array<OneD, NekDouble> wsp(nquad1 * order0);
        IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
                                     m_base[1]->GetBdata(), inarray, outarray,
                                     wsp);
    }
}
 
void StdTriExp::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)
{
    boost::ignore_unused(doCheckCollDir0, doCheckCollDir1);
 
    int i;
    int mode;
    int nquad0  = m_base[0]->GetNumPoints();
    int nquad1  = m_base[1]->GetNumPoints();
    int nmodes0 = m_base[0]->GetNumModes();
    int nmodes1 = m_base[1]->GetNumModes();
 
    ASSERTL1(wsp.size() >= nquad1 * nmodes0,
             "Workspace size is not sufficient");
 
    Blas::Dgemm('T', 'N', nquad1, nmodes0, nquad0, 1.0, inarray.get(), nquad0,
                base0.get(), nquad0, 0.0, wsp.get(), nquad1);
 
    // Inner product with respect to 'b' direction
    for (mode = i = 0; i < nmodes0; ++i)
    {
        Blas::Dgemv('T', nquad1, nmodes1 - i, 1.0, base1.get() + mode * nquad1,
                    nquad1, wsp.get() + i * nquad1, 1, 0.0,
                    outarray.get() + mode, 1);
        mode += nmodes1 - i;
    }
 
    // fix for modified basis by splitting top vertex mode
    if (m_base[0]->GetBasisType() == LibUtilities::eModified_A)
    {
        outarray[1] +=
            Blas::Ddot(nquad1, base1.get() + nquad1, 1, wsp.get() + nquad1, 1);
    }
}
 
void StdTriExp::v_IProductWRTDerivBase(
    const int dir, const Array<OneD, const NekDouble> &inarray,
    Array<OneD, NekDouble> &outarray)
{
    StdTriExp::v_IProductWRTDerivBase_SumFac(dir, inarray, outarray);
}
 
void StdTriExp::v_IProductWRTDerivBase_SumFac(
    const int dir, const Array<OneD, const NekDouble> &inarray,
    Array<OneD, NekDouble> &outarray)
{
    int i;
    int nquad0  = m_base[0]->GetNumPoints();
    int nquad1  = m_base[1]->GetNumPoints();
    int nqtot   = nquad0 * nquad1;
    int nmodes0 = m_base[0]->GetNumModes();
    int wspsize = max(max(nqtot, m_ncoeffs), nquad1 * nmodes0);
 
    Array<OneD, NekDouble> gfac0(2 * wspsize);
    Array<OneD, NekDouble> tmp0(gfac0 + wspsize);
 
    const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
 
    // set up geometric factor: 2/(1-z1)
    for (i = 0; i < nquad1; ++i)
    {
        gfac0[i] = 2.0 / (1 - z1[i]);
    }
 
    for (i = 0; i < nquad1; ++i)
    {
        Vmath::Smul(nquad0, gfac0[i], &inarray[0] + i * nquad0, 1,
                    &tmp0[0] + i * nquad0, 1);
    }
 
    MultiplyByQuadratureMetric(tmp0, tmp0);
 
    switch (dir)
    {
        case 0:
        {
            IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
                                         m_base[1]->GetBdata(), tmp0, outarray,
                                         gfac0);
            break;
        }
        case 1:
        {
            Array<OneD, NekDouble> tmp3(m_ncoeffs);
            const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
 
            for (i = 0; i < nquad0; ++i)
            {
                gfac0[i] = 0.5 * (1 + z0[i]);
            }
 
            for (i = 0; i < nquad1; ++i)
            {
                Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp0[0] + i * nquad0, 1,
                            &tmp0[0] + i * nquad0, 1);
            }
 
            IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
                                         m_base[1]->GetBdata(), tmp0, tmp3,
                                         gfac0);
 
            MultiplyByQuadratureMetric(inarray, tmp0);
            IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
                                         m_base[1]->GetDbdata(), tmp0, outarray,
                                         gfac0);
            Vmath::Vadd(m_ncoeffs, &tmp3[0], 1, &outarray[0], 1, &outarray[0],
                        1);
            break;
        }
        default:
        {
            ASSERTL1(false, "input dir is out of range");
            break;
        }
    }
}
 
//---------------------------------------
// Evaluation functions
//---------------------------------------
 
void StdTriExp::v_LocCoordToLocCollapsed(const Array<OneD, const NekDouble> &xi,
                                         Array<OneD, NekDouble> &eta)
{
    NekDouble d1 = 1. - xi[1];
    if (fabs(d1) < NekConstants::kNekZeroTol)
    {
        if (d1 >= 0.)
        {
            d1 = NekConstants::kNekZeroTol;
        }
        else
        {
            d1 = -NekConstants::kNekZeroTol;
        }
    }
    eta[0] = 2. * (1. + xi[0]) / d1 - 1.0;
    eta[1] = xi[1];
}
 
void StdTriExp::v_LocCollapsedToLocCoord(
    const Array<OneD, const NekDouble> &eta, Array<OneD, NekDouble> &xi)
{
    xi[0] = (1.0 + eta[0]) * (1.0 - eta[1]) * 0.5 - 1.0;
    xi[1] = eta[1];
}
 
void StdTriExp::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
{
    int i, m;
    int nquad0                         = m_base[0]->GetNumPoints();
    int nquad1                         = m_base[1]->GetNumPoints();
    int order0                         = m_base[0]->GetNumModes();
    int order1                         = m_base[1]->GetNumModes();
    int mode0                          = 0;
    Array<OneD, const NekDouble> base0 = m_base[0]->GetBdata();
    Array<OneD, const NekDouble> base1 = m_base[1]->GetBdata();
 
    ASSERTL2(mode <= m_ncoeffs, "calling argument mode is larger than "
                                "total expansion order");
 
    m = order1;
    for (i = 0; i < order0; ++i, m += order1 - i)
    {
        if (m > mode)
        {
            mode0 = i;
            break;
        }
    }
 
    // deal with top vertex mode in modified basis
    if (mode == 1 && m_base[0]->GetBasisType() == LibUtilities::eModified_A)
    {
        Vmath::Fill(nquad0 * nquad1, 1.0, outarray, 1);
    }
    else
    {
        for (i = 0; i < nquad1; ++i)
        {
            Vmath::Vcopy(nquad0, (NekDouble *)(base0.get() + mode0 * nquad0), 1,
                         &outarray[0] + i * nquad0, 1);
        }
    }
 
    for (i = 0; i < nquad0; ++i)
    {
        Vmath::Vmul(nquad1, (NekDouble *)(base1.get() + mode * nquad1), 1,
                    &outarray[0] + i, nquad0, &outarray[0] + i, nquad0);
    }
}
 
NekDouble StdTriExp::v_PhysEvaluateBasis(
    const Array<OneD, const NekDouble> &coords, int mode)
{
    Array<OneD, NekDouble> coll(2);
    LocCoordToLocCollapsed(coords, coll);
 
    // From mode we need to determine mode0 and mode1 in the (p,q)
    // direction. mode1 can be directly inferred from mode.
    const int nm1   = m_base[1]->GetNumModes();
    const double c  = 1 + 2 * nm1;
    const int mode0 = floor(0.5 * (c - sqrt(c * c - 8 * mode)));
 
    if (mode == 1 && m_base[0]->GetBasisType() == LibUtilities::eModified_A)
    {
        // Account for collapsed vertex.
        return StdExpansion::BaryEvaluateBasis<1>(coll[1], 1);
    }
    else
    {
        return StdExpansion::BaryEvaluateBasis<0>(coll[0], mode0) *
               StdExpansion::BaryEvaluateBasis<1>(coll[1], mode);
    }
}
 
NekDouble StdTriExp::v_PhysEvaluate(const Array<OneD, NekDouble> &coord,
                                    const Array<OneD, const NekDouble> &inarray,
                                    std::array<NekDouble, 3> &firstOrderDerivs)
{
    // Collapse coordinates
    Array<OneD, NekDouble> coll(2, 0.0);
    LocCoordToLocCollapsed(coord, coll);
 
    // If near singularity do the old interpolation matrix method
    if ((1 - coll[1]) < 1e-5)
    {
        int totPoints = GetTotPoints();
        Array<OneD, NekDouble> EphysDeriv0(totPoints), EphysDeriv1(totPoints);
        PhysDeriv(inarray, EphysDeriv0, EphysDeriv1);
 
        Array<OneD, DNekMatSharedPtr> I(2);
        I[0] = GetBase()[0]->GetI(coll);
        I[1] = GetBase()[1]->GetI(coll + 1);
 
        firstOrderDerivs[0] = PhysEvaluate(I, EphysDeriv0);
        firstOrderDerivs[1] = PhysEvaluate(I, EphysDeriv1);
        return PhysEvaluate(I, inarray);
    }
 
    // set up geometric factor: 2.0/(1.0-z1)
    NekDouble fac0 = 2 / (1 - coll[1]);
 
    NekDouble val = BaryTensorDeriv(coll, inarray, firstOrderDerivs);
 
    // Copy d0 into temp for d1
    NekDouble temp;
    temp = firstOrderDerivs[0];
 
    // Multiply by geometric factor
    firstOrderDerivs[0] = firstOrderDerivs[0] * fac0;
 
    // set up geometric factor: (1+z0)/(1-z1)
    NekDouble fac1 = fac0 * (coll[0] + 1) / 2;
 
    // Multiply out_d0 by geometric factor and add to out_d1
    firstOrderDerivs[1] += fac1 * temp;
 
    return val;
}
 
int StdTriExp::v_GetNverts() const
{
    return 3;
}
 
int StdTriExp::v_GetNtraces() const
{
    return 3;
}
 
LibUtilities::ShapeType StdTriExp::v_DetShapeType() const
{
    return LibUtilities::eTriangle;
}
 
int StdTriExp::v_NumBndryCoeffs() const
{
    ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A,
             "BasisType is not a boundary interior form");
    ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B,
             "BasisType is not a boundary interior form");
 
    return 3 + (GetBasisNumModes(0) - 2) + 2 * (GetBasisNumModes(1) - 2);
}
 
int StdTriExp::v_NumDGBndryCoeffs() const
{
    ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A,
             "BasisType is not a boundary interior form");
    ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B,
             "BasisType is not a boundary interior form");
 
    return GetBasisNumModes(0) + 2 * GetBasisNumModes(1);
}
 
int StdTriExp::v_GetTraceNcoeffs(const int i) const
{
    ASSERTL2(i >= 0 && i <= 2, "edge id is out of range");
 
    if (i == 0)
    {
        return GetBasisNumModes(0);
    }
    else
    {
        return GetBasisNumModes(1);
    }
}
 
int StdTriExp::v_GetTraceIntNcoeffs(const int i) const
{
    ASSERTL2(i >= 0 && i <= 2, "edge id is out of range");
 
    if (i == 0)
    {
        return GetBasisNumModes(0) - 2;
    }
    else
    {
        return GetBasisNumModes(1) - 2;
    }
}
 
int StdTriExp::v_GetTraceNumPoints(const int i) const
{
    ASSERTL2((i >= 0) && (i <= 2), "edge id is out of range");
 
    if (i == 0)
    {
        return GetNumPoints(0);
    }
    else
    {
        return GetNumPoints(1);
    }
}
 
int StdTriExp::v_CalcNumberOfCoefficients(
    const std::vector<unsigned int> &nummodes, int &modes_offset)
{
    int nmodes = LibUtilities::StdTriData::getNumberOfCoefficients(
        nummodes[modes_offset], nummodes[modes_offset + 1]);
    modes_offset += 2;
 
    return nmodes;
}
 
void StdTriExp::v_GetCoords(Array<OneD, NekDouble> &coords_0,
                            Array<OneD, NekDouble> &coords_1,
                            Array<OneD, NekDouble> &coords_2)
{
    boost::ignore_unused(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, j;
 
    for (i = 0; i < nq1; ++i)
    {
        for (j = 0; j < nq0; ++j)
        {
            coords_0[i * nq0 + j] = (1 + z0[j]) * (1 - z1[i]) / 2.0 - 1.0;
        }
        Vmath::Fill(nq0, z1[i], &coords_1[0] + i * nq0, 1);
    }
}
 
bool StdTriExp::v_IsBoundaryInteriorExpansion() const
{
    return m_base[0]->GetBasisType() == LibUtilities::eModified_A &&
           m_base[1]->GetBasisType() == LibUtilities::eModified_B;
}
 
const LibUtilities::BasisKey StdTriExp::v_GetTraceBasisKey(const int i,
                                                           const int j) const
{
    boost::ignore_unused(j);
    ASSERTL2(i >= 0 && i <= 2, "edge id is out of range");
 
    if (i == 0)
    {
        return GetBasis(0)->GetBasisKey();
    }
    else
    {
        switch (m_base[1]->GetBasisType())
        {
            case LibUtilities::eModified_B:
            {
                switch (m_base[1]->GetPointsType())
                {
                    case LibUtilities::eGaussRadauMAlpha1Beta0:
                    {
                        LibUtilities::PointsKey pkey(
                            m_base[1]
                                    ->GetBasisKey()
                                    .GetPointsKey()
                                    .GetNumPoints() +
                                1,
                            LibUtilities::eGaussLobattoLegendre);
                        return LibUtilities::BasisKey(LibUtilities::eModified_A,
                                                      m_base[1]->GetNumModes(),
                                                      pkey);
                        break;
                    }
 
                    break;
                        // Currently this does not increase the points by
                        // 1 since when using this quadrature we are
                        // presuming it is already been increased by one
                        // when comopared to
                        // GaussRadauMAlpha1Beta0. Currently used in the
                        // GJP option
                        //
                        // Note have put down it back to numpoints +1 to
                        // test for use on tri faces and GJP.
 
                    case LibUtilities::eGaussRadauMLegendre:
                    {
                        LibUtilities::PointsKey pkey(
                            m_base[1]
                                    ->GetBasisKey()
                                    .GetPointsKey()
                                    .GetNumPoints() +
                                1,
                            LibUtilities::eGaussLobattoLegendre);
                        return LibUtilities::BasisKey(LibUtilities::eModified_A,
                                                      m_base[1]->GetNumModes(),
                                                      pkey);
                        break;
                    }
 
                    case LibUtilities::ePolyEvenlySpaced:
                    {
                        LibUtilities::PointsKey pkey(
                            m_base[1]
                                    ->GetBasisKey()
                                    .GetPointsKey()
                                    .GetNumPoints() +
                                1,
                            LibUtilities::ePolyEvenlySpaced);
                        return LibUtilities::BasisKey(LibUtilities::eModified_A,
                                                      m_base[1]->GetNumModes(),
                                                      pkey);
                        break;
                    }
 
                    break;
                    default:
                        NEKERROR(ErrorUtil::efatal,
                                 "Unexpected points distribution " +
                                     LibUtilities::kPointsTypeStr
                                         [m_base[1]->GetPointsType()] +
                                     " in StdTriExp::v_GetTraceBasisKey");
                        break;
                }
                break;
            }
            default:
                NEKERROR(ErrorUtil::efatal,
                         "Information not available to set edge key");
                break;
        }
    }
    return LibUtilities::NullBasisKey;
}
 
//--------------------------
// Mappings
//--------------------------
 
int StdTriExp::v_GetVertexMap(const int localVertexId, bool useCoeffPacking)
{
    ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
                 GetBasisType(1) == LibUtilities::eModified_B,
             "Mapping not defined for this type of basis");
 
    int localDOF = 0;
    if (useCoeffPacking == true)
    {
        switch (localVertexId)
        {
            case 0:
            {
                localDOF = 0;
                break;
            }
            case 1:
            {
                localDOF = 1;
                break;
            }
            case 2:
            {
                localDOF = m_base[1]->GetNumModes();
                break;
            }
            default:
            {
                ASSERTL0(false, "eid must be between 0 and 2");
                break;
            }
        }
    }
    else // follow book format for vertex indexing.
    {
        switch (localVertexId)
        {
            case 0:
            {
                localDOF = 0;
                break;
            }
            case 1:
            {
                localDOF = m_base[1]->GetNumModes();
                break;
            }
            case 2:
            {
                localDOF = 1;
                break;
            }
            default:
            {
                ASSERTL0(false, "eid must be between 0 and 2");
                break;
            }
        }
    }
 
    return localDOF;
}
 
void StdTriExp::v_GetInteriorMap(Array<OneD, unsigned int> &outarray)
{
    ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A &&
                 GetBasisType(1) == LibUtilities::eModified_B,
             "Expansion not of a proper type");
 
    int i, j;
    int cnt = 0;
    int nummodes0, nummodes1;
    int startvalue;
    if (outarray.size() != GetNcoeffs() - NumBndryCoeffs())
    {
        outarray = Array<OneD, unsigned int>(GetNcoeffs() - NumBndryCoeffs());
    }
 
    nummodes0 = m_base[0]->GetNumModes();
    nummodes1 = m_base[1]->GetNumModes();
 
    startvalue = 2 * nummodes1;
 
    for (i = 0; i < nummodes0 - 2; i++)
    {
        for (j = 0; j < nummodes1 - 3 - i; j++)
        {
            outarray[cnt++] = startvalue + j;
        }
        startvalue += nummodes1 - 2 - i;
    }
}
 
void StdTriExp::v_GetBoundaryMap(Array<OneD, unsigned int> &outarray)
{
    ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A &&
                 GetBasisType(1) == LibUtilities::eModified_B,
             "Expansion not of expected type");
    int i;
    int cnt;
    int nummodes0, nummodes1;
    int value;
 
    if (outarray.size() != NumBndryCoeffs())
    {
        outarray = Array<OneD, unsigned int>(NumBndryCoeffs());
    }
 
    nummodes0 = m_base[0]->GetNumModes();
    nummodes1 = m_base[1]->GetNumModes();
 
    value = 2 * nummodes1 - 1;
    for (i = 0; i < value; i++)
    {
        outarray[i] = i;
    }
    cnt = value;
 
    for (i = 0; i < nummodes0 - 2; i++)
    {
        outarray[cnt++] = value;
        value += nummodes1 - 2 - i;
    }
}
 
void StdTriExp::v_GetTraceCoeffMap(const unsigned int eid,
                                   Array<OneD, unsigned int> &maparray)
{
 
    ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
                 GetBasisType(1) == LibUtilities::eModified_B,
             "Mapping not defined for this type of basis");
 
    ASSERTL1(eid < 3, "eid must be between 0 and 2");
 
    int i;
    int order0   = m_base[0]->GetNumModes();
    int order1   = m_base[1]->GetNumModes();
    int numModes = (eid == 0) ? order0 : order1;
 
    if (maparray.size() != numModes)
    {
        maparray = Array<OneD, unsigned int>(numModes);
    }
 
    switch (eid)
    {
        case 0:
        {
            int cnt = 0;
            for (i = 0; i < numModes; cnt += order1 - i, ++i)
            {
                maparray[i] = cnt;
            }
            break;
        }
        case 1:
        {
            maparray[0] = order1;
            maparray[1] = 1;
            for (i = 2; i < numModes; i++)
            {
                maparray[i] = order1 - 1 + i;
            }
            break;
        }
        case 2:
        {
            for (i = 0; i < numModes; i++)
            {
                maparray[i] = i;
            }
            break;
        }
        default:
            ASSERTL0(false, "eid must be between 0 and 2");
            break;
    }
}
 
void StdTriExp::v_GetTraceInteriorToElementMap(
    const int eid, Array<OneD, unsigned int> &maparray,
    Array<OneD, int> &signarray, const Orientation edgeOrient)
{
    ASSERTL0(GetBasisType(0) == LibUtilities::eModified_A ||
                 GetBasisType(1) == LibUtilities::eModified_B,
             "Mapping not defined for this type of basis");
    int i;
    const int nummodes1      = m_base[1]->GetNumModes();
    const int nEdgeIntCoeffs = GetTraceNcoeffs(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.get(), signarray.get() + nEdgeIntCoeffs, 1);
    }
 
    switch (eid)
    {
        case 0:
        {
            int cnt = 2 * nummodes1 - 1;
            for (i = 0; i < nEdgeIntCoeffs; cnt += nummodes1 - 2 - i, ++i)
            {
                maparray[i] = cnt;
            }
            break;
        }
        case 1:
        {
            for (i = 0; i < nEdgeIntCoeffs; i++)
            {
                maparray[i] = nummodes1 + 1 + i;
            }
            break;
        }
        case 2:
        {
            for (i = 0; i < nEdgeIntCoeffs; i++)
            {
                maparray[i] = 2 + i;
            }
            break;
        }
        default:
        {
            ASSERTL0(false, "eid must be between 0 and 2");
            break;
        }
    }
 
    if (edgeOrient == eBackwards)
    {
        for (i = 1; i < nEdgeIntCoeffs; i += 2)
        {
            signarray[i] = -1;
        }
    }
}
 
//---------------------------------------
// Wrapper functions
//---------------------------------------
 
DNekMatSharedPtr StdTriExp::v_GenMatrix(const StdMatrixKey &mkey)
{
 
    MatrixType mtype = mkey.GetMatrixType();
 
    DNekMatSharedPtr Mat;
 
    switch (mtype)
    {
        case ePhysInterpToEquiSpaced:
        {
            int nq0, nq1, nq;
 
            nq0 = m_base[0]->GetNumPoints();
            nq1 = m_base[1]->GetNumPoints();
 
            // take definition from key
            if (mkey.ConstFactorExists(eFactorConst))
            {
                nq = (int)mkey.GetConstFactor(eFactorConst);
            }
            else
            {
                nq = max(nq0, nq1);
            }
 
            int neq = LibUtilities::StdTriData::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 (int i = 0; i < nq; ++i)
            {
                for (int j = 0; j < nq - i; ++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 (int 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 (int 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;
        }
        default:
        {
            Mat = StdExpansion::CreateGeneralMatrix(mkey);
            break;
        }
    }
 
    return Mat;
}
 
DNekMatSharedPtr StdTriExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
{
    return v_GenMatrix(mkey);
}
 
//---------------------------------------
// Operator evaluation functions
//---------------------------------------
 
void StdTriExp::v_MassMatrixOp(const Array<OneD, const NekDouble> &inarray,
                               Array<OneD, NekDouble> &outarray,
                               const StdMatrixKey &mkey)
{
    StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
}
 
void StdTriExp::v_LaplacianMatrixOp(const Array<OneD, const NekDouble> &inarray,
                                    Array<OneD, NekDouble> &outarray,
                                    const StdMatrixKey &mkey)
{
    StdTriExp::v_LaplacianMatrixOp_MatFree(inarray, outarray, mkey);
}
 
void StdTriExp::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 StdTriExp::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 StdTriExp::v_HelmholtzMatrixOp(const Array<OneD, const NekDouble> &inarray,
                                    Array<OneD, NekDouble> &outarray,
                                    const StdMatrixKey &mkey)
{
    StdTriExp::v_HelmholtzMatrixOp_MatFree(inarray, outarray, mkey);
}
 
void StdTriExp::v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
                                     const StdMatrixKey &mkey)
{
    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();
 
    // 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_B, nmodes_b, pb);
    StdTriExp OrthoExp(Ba, Bb);
 
    Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
 
    // project onto physical space.
    OrthoExp.FwdTrans(array, orthocoeffs);
 
    if (mkey.ConstFactorExists(
            eFactorSVVPowerKerDiffCoeff)) // Rodrigo's power kern
    {
        NekDouble cutoff = mkey.GetConstFactor(eFactorSVVCutoffRatio);
        NekDouble SvvDiffCoeff =
            mkey.GetConstFactor(eFactorSVVPowerKerDiffCoeff) *
            mkey.GetConstFactor(eFactorSVVDiffCoeff);
 
        int cnt = 0;
        for (int j = 0; j < nmodes_a; ++j)
        {
            for (int k = 0; k < nmodes_b - j; ++k, ++cnt)
            {
                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[cnt] *= (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);
 
        int cnt = 0;
        for (int j = 0; j < nmodes_a; ++j)
        {
            for (int k = 0; k < nmodes_b - j; ++k, ++cnt)
            {
                int maxjk = max(j, k);
                maxjk     = min(maxjk, kSVVDGFiltermodesmax - 1);
 
                orthocoeffs[cnt] *= SvvDiffCoeff * kSVVDGFilter[max_ab][maxjk];
            }
        }
    }
    else
    {
        NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDiffCoeff);
 
        int cutoff = (int)(mkey.GetConstFactor(eFactorSVVCutoffRatio) *
                           min(nmodes_a, nmodes_b));
 
        NekDouble epsilon = 1.0;
        int nmodes        = min(nmodes_a, nmodes_b);
 
        int cnt = 0;
 
        // apply SVV filter (JEL)
        for (int j = 0; j < nmodes_a; ++j)
        {
            for (int k = 0; k < nmodes_b - j; ++k)
            {
                if (j + k >= cutoff)
                {
                    orthocoeffs[cnt] *=
                        (SvvDiffCoeff *
                         exp(-(j + k - nmodes) * (j + k - nmodes) /
                             ((NekDouble)((j + k - cutoff + epsilon) *
                                          (j + k - cutoff + epsilon)))));
                }
                else
                {
                    orthocoeffs[cnt] *= 0.0;
                }
                cnt++;
            }
        }
    }
 
    // backward transform to physical space
    OrthoExp.BwdTrans(orthocoeffs, array);
}
 
void StdTriExp::v_ReduceOrderCoeffs(int numMin,
                                    const Array<OneD, const NekDouble> &inarray,
                                    Array<OneD, NekDouble> &outarray)
{
    int n_coeffs = inarray.size();
    int nquad0   = m_base[0]->GetNumPoints();
    int nquad1   = m_base[1]->GetNumPoints();
    Array<OneD, NekDouble> coeff(n_coeffs);
    Array<OneD, NekDouble> coeff_tmp(n_coeffs, 0.0);
    Array<OneD, NekDouble> tmp;
    Array<OneD, NekDouble> tmp2;
    int nqtot = nquad0 * nquad1;
    Array<OneD, NekDouble> phys_tmp(nqtot, 0.0);
 
    int nmodes0 = m_base[0]->GetNumModes();
    int nmodes1 = m_base[1]->GetNumModes();
    int numMin2 = nmodes0;
    int i;
 
    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_B, nmodes1, Pkey1);
 
    StdRegions::StdTriExpSharedPtr m_OrthoTriExp;
    StdRegions::StdTriExpSharedPtr m_TriExp;
 
    m_TriExp = MemoryManager<StdRegions::StdTriExp>::AllocateSharedPtr(b0, b1);
    m_OrthoTriExp = MemoryManager<StdRegions::StdTriExp>::AllocateSharedPtr(
        bortho0, bortho1);
 
    m_TriExp->BwdTrans(inarray, phys_tmp);
    m_OrthoTriExp->FwdTrans(phys_tmp, coeff);
 
    for (i = 0; i < n_coeffs; i++)
    {
        if (i == numMin)
        {
            coeff[i] = 0.0;
            numMin += numMin2 - 1;
            numMin2 -= 1.0;
        }
    }
 
    m_OrthoTriExp->BwdTrans(coeff, phys_tmp);
    m_TriExp->FwdTrans(phys_tmp, outarray);
}
 
//---------------------------------------
// Private helper functions
//---------------------------------------
 
void StdTriExp::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();
    const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
 
    // multiply by integration constants
    for (i = 0; i < nquad1; ++i)
    {
        Vmath::Vmul(nquad0, inarray.get() + i * nquad0, 1, w0.get(), 1,
                    outarray.get() + i * nquad0, 1);
    }
 
    switch (m_base[1]->GetPointsType())
    {
            // (1,0) Jacobi Inner product
        case LibUtilities::eGaussRadauMAlpha1Beta0:
            for (i = 0; i < nquad1; ++i)
            {
                Blas::Dscal(nquad0, 0.5 * w1[i], outarray.get() + i * nquad0,
                            1);
            }
            break;
            // Legendre inner product
        default:
            for (i = 0; i < nquad1; ++i)
            {
                Blas::Dscal(nquad0, 0.5 * (1 - z1[i]) * w1[i],
                            outarray.get() + i * nquad0, 1);
            }
            break;
    }
}
 
void StdTriExp::v_GetSimplexEquiSpacedConnectivity(Array<OneD, int> &conn,
                                                   bool standard)
{
    boost::ignore_unused(standard);
 
    int np1 = m_base[0]->GetNumPoints();
    int np2 = m_base[1]->GetNumPoints();
    int np  = max(np1, np2);
 
    conn = Array<OneD, int>(3 * (np - 1) * (np - 1));
 
    int row   = 0;
    int rowp1 = 0;
    int cnt   = 0;
    for (int i = 0; i < np - 1; ++i)
    {
        rowp1 += np - i;
        for (int j = 0; j < np - i - 2; ++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;
        }
 
        conn[cnt++] = row + np - i - 2;
        conn[cnt++] = row + np - i - 1;
        conn[cnt++] = rowp1 + np - i - 2;
 
        row += np - i;
    }
}
} // namespace StdRegions
} // namespace Nektar