StdExpansion.h

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File: StdExpansion.h
//
// 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: Class definition StdExpansion which is the base class
// to all expansion shapes
//
///////////////////////////////////////////////////////////////////////////////
 
#ifndef NEKTAR_LIB_STDREGIONS_STANDARDEXPANSION_H
#define NEKTAR_LIB_STDREGIONS_STANDARDEXPANSION_H
 
#include <fstream>
#include <memory>
#include <vector>
 
#include <boost/core/ignore_unused.hpp>
#include <set>
 
#include <LibUtilities/LinearAlgebra/NekTypeDefs.hpp>
#include <StdRegions/StdMatrixKey.h>
#include <StdRegions/StdRegions.hpp>
#include <StdRegions/StdRegionsDeclspec.h>
namespace Nektar
{
namespace LocalRegions
{
class MatrixKey;
class Expansion;
} // namespace LocalRegions
} // namespace Nektar
 
namespace Nektar
{
namespace StdRegions
{
 
/** \brief The base class for all shapes
 *
 *  This is the lowest level basic class for all shapes and so
 *  contains the definition of common data and common routine to all
 *  elements
 */
class StdExpansion : public std::enable_shared_from_this<StdExpansion>
{
public:
    /** \brief Default Constructor */
    STD_REGIONS_EXPORT StdExpansion();
 
    /** \brief Constructor */
    STD_REGIONS_EXPORT StdExpansion(
        const int numcoeffs, const int numbases,
        const LibUtilities::BasisKey &Ba = LibUtilities::NullBasisKey,
        const LibUtilities::BasisKey &Bb = LibUtilities::NullBasisKey,
        const LibUtilities::BasisKey &Bc = LibUtilities::NullBasisKey);
 
    /** \brief Copy Constructor */
    STD_REGIONS_EXPORT StdExpansion(const StdExpansion &T);
 
    /** \brief Destructor */
    STD_REGIONS_EXPORT virtual ~StdExpansion();
 
    // Standard Expansion Routines Applicable Regardless of Region
 
    /** \brief This function returns the number of 1D bases used in
     *  the expansion
     *
     *  \return returns the number of 1D bases used in the expansion,
     *  which is equal to number dimension of the expansion
     */
    inline int GetNumBases() const
    {
        return m_base.size();
    }
 
    /** \brief This function gets the shared point to basis
     *
     *  \return returns the shared pointer to the bases
     */
    inline const Array<OneD, const LibUtilities::BasisSharedPtr> &GetBase()
        const
    {
        return (m_base);
    }
 
    /** \brief This function gets the shared point to basis in
     *  the \a dir direction
     *
     *  \return returns the shared pointer to the basis in
     *  directin \a dir
     */
    inline const LibUtilities::BasisSharedPtr &GetBasis(int dir) const
    {
        ASSERTL1(dir < m_base.size(), "dir is larger than number of bases");
        return (m_base[dir]);
    }
 
    /** \brief This function returns the total number of coefficients
     *  used in the expansion
     *
     *  \return returns the total number of coefficients (which is
     *  equivalent to the total number of modes) used in the expansion
     */
    inline int GetNcoeffs(void) const
    {
        return (m_ncoeffs);
    }
 
    /** \brief This function returns the total number of quadrature
     *  points used in the element
     *
     *  \return returns the total number of quadrature points
     */
    inline int GetTotPoints() const
    {
        int nqtot = 1;
 
        for (size_t i = 0; i < m_base.size(); ++i)
        {
            nqtot *= m_base[i]->GetNumPoints();
        }
 
        return nqtot;
    }
 
    /** \brief This function returns the type of basis used in the \a dir
     *  direction
     *
     *  The different types of bases implemented in the code are defined
     *  in the LibUtilities::BasisType enumeration list. As a result, the
     *  function will return one of the types of this enumeration list.
     *
     *  \param dir the direction
     *  \return returns the type of basis used in the \a dir direction
     */
    inline LibUtilities::BasisType GetBasisType(const int dir) const
    {
        ASSERTL1(dir < m_base.size(), "dir is larger than m_numbases");
        return (m_base[dir]->GetBasisType());
    }
 
    /** \brief This function returns the number of expansion modes
     *  in the \a dir direction
     *
     *  \param dir the direction
     *  \return returns the number of expansion modes in the \a dir
     *  direction
     */
    inline int GetBasisNumModes(const int dir) const
    {
        ASSERTL1(dir < m_base.size(), "dir is larger than m_numbases");
        return (m_base[dir]->GetNumModes());
    }
 
    /** \brief This function returns the maximum number of
     *  expansion modes over all local directions
     *
     *  \return returns the maximum number of expansion modes
     *  over all local directions
     */
    inline int EvalBasisNumModesMax(void) const
    {
        int returnval = 0;
 
        for (size_t i = 0; i < m_base.size(); ++i)
        {
            returnval = std::max(returnval, m_base[i]->GetNumModes());
        }
 
        return returnval;
    }
 
    /** \brief This function returns the type of quadrature points used
     *  in the \a dir direction
     *
     *  The different types of quadrature points implemented in the code
     *  are defined in the LibUtilities::PointsType enumeration list.
     *  As a result, the function will return one of the types of this
     *  enumeration list.
     *
     *  \param dir the direction
     *  \return returns the type of quadrature points  used in the \a dir
     *  direction
     */
    inline LibUtilities::PointsType GetPointsType(const int dir) const
    {
        ASSERTL1(dir < m_base.size(), "dir is larger than m_numbases");
        return (m_base[dir]->GetPointsType());
    }
 
    /** \brief This function returns the number of quadrature points
     *  in the \a dir direction
     *
     *  \param dir the direction
     *  \return returns the number of quadrature points in the \a dir
     *  direction
     */
    inline int GetNumPoints(const int dir) const
    {
        ASSERTL1(dir < m_base.size() || dir == 0,
                 "dir is larger than m_numbases");
        return (m_base.size() > 0 ? m_base[dir]->GetNumPoints() : 1);
    }
 
    /** \brief This function returns a pointer to the array containing
     *  the quadrature points in \a dir direction
     *
     *  \param dir the direction
     *  \return returns a pointer to the array containing
     *  the quadrature points in \a dir direction
     */
    inline const Array<OneD, const NekDouble> &GetPoints(const int dir) const
    {
        return m_base[dir]->GetZ();
    }
 
    // Wrappers around virtual Functions
    /** \brief This function returns the number of vertices of the
     *  expansion domain
     *
     *  This function is a wrapper around the virtual function
     *  \a v_GetNverts()
     *
     *  \return returns the number of vertices of the expansion domain
     */
    int GetNverts() const
    {
        return v_GetNverts();
    }
 
    /** \brief This function returns the number of expansion coefficients
     *  belonging to the \a i-th trace
     *
     *  This function is a wrapper around the virtual function
     *  \a v_GetTraceNcoeffs()
     *
     *  \param i specifies which trace
     *  \return returns the number of expansion coefficients belonging to
     *  the \a i-th trace
     */
    int GetTraceNcoeffs(const int i) const
    {
        return v_GetTraceNcoeffs(i);
    }
 
    int GetTraceIntNcoeffs(const int i) const
    {
        return v_GetTraceIntNcoeffs(i);
    }
 
    /** \brief This function returns the number of quadrature points
     *  belonging to the \a i-th trace
     *
     *  This function is a wrapper around the virtual function
     *  \a v_GetTraceNumPoints()
     *
     *  \param i specifies which trace id
     *  \return returns the number of quadrature points belonging to
     *   the \a i-th trace
     */
    int GetTraceNumPoints(const int i) const
    {
        return v_GetTraceNumPoints(i);
    }
 
    /** \brief This function returns the basis key belonging
     *   to the \a i-th trace
     *
     *  This function is a wrapper around the virtual function
     *  \a v_GetTraceBasisKey()
     *
     *  \param i specifies which trace id
     *  \param k is the direction of the basis key for 2D traces
     *
     *  \return returns the number of Basis key of the ith
     *  trace in the k th direction (when trace is a 2D
     *  object)
     */
    const LibUtilities::BasisKey GetTraceBasisKey(const int i, int k = -1) const
    {
        return v_GetTraceBasisKey(i, k);
    }
 
    /** \brief This function returns the basis key belonging
     *   to the \a i-th trace
     *
     *  This function is a wrapper around the virtual function
     *  \a v_GetTracePointsKey()
     *
     *  \param i specifies which trace id
     *  \param k is the direction of the basis key for 2D traces
     *
     *  \return returns the number of Points key of the ith
     *  trace in the k th direction (when trace is a 2D
     *  object)
     */
    LibUtilities::PointsKey GetTracePointsKey(const int i, int k = -1) const
    {
        return v_GetTracePointsKey(i, k);
    }
 
    int NumBndryCoeffs(void) const
    {
        return v_NumBndryCoeffs();
    }
 
    int NumDGBndryCoeffs(void) const
    {
        return v_NumDGBndryCoeffs();
    }
 
    /** \brief This function returns the type of expansion
     *  Nodal point type if defined
     *
     *  This function is a wrapper around the virtual function
     *  \a v_GetNodalPointsKey()
     *
     */
    const LibUtilities::PointsKey GetNodalPointsKey() const
    {
        return v_GetNodalPointsKey();
    };
 
    /**
     * @brief Returns the number of trace elements connected to this
     * element.
     *
     * For example, a quadrilateral has four edges, so this function
     * would return 4.
     */
    int GetNtraces() const
    {
        return v_GetNtraces();
    }
 
    /** \brief This function returns the shape of the expansion domain
     *
     *  This function is a wrapper around the virtual function
     *  \a v_DetShapeType()
     *
     *  The different shape types implemented in the code are defined
     *  in the ::ShapeType enumeration list. As a result, the
     *  function will return one of the types of this enumeration list.
     *
     *  \return returns the shape of the expansion domain
     */
    LibUtilities::ShapeType DetShapeType() const
    {
        return v_DetShapeType();
    }
 
    std::shared_ptr<StdExpansion> GetStdExp() const
    {
        return v_GetStdExp();
    }
 
    std::shared_ptr<StdExpansion> GetLinStdExp(void) const
    {
        return v_GetLinStdExp();
    }
 
    int GetShapeDimension() const
    {
        return v_GetShapeDimension();
    }
 
    bool IsBoundaryInteriorExpansion() const
    {
        return v_IsBoundaryInteriorExpansion();
    }
 
    bool IsNodalNonTensorialExp()
    {
        return v_IsNodalNonTensorialExp();
    }
 
    /** \brief This function performs the Backward transformation from
     *  coefficient space to physical space
     *
     *  This function is a wrapper around the virtual function
     *  \a v_BwdTrans()
     *
     *  Based on the expansion coefficients, this function evaluates the
     *  expansion at the quadrature points. This is equivalent to the
     *  operation \f[ u(\xi_{1i}) =
     *  \sum_{p=0}^{P-1} \hat{u}_p \phi_p(\xi_{1i}) \f] which can be
     *  evaluated as \f$ {\bf u} = {\bf B}^T {\bf \hat{u}} \f$ with
     *  \f${\bf B}[i][j] = \phi_i(\xi_{j})\f$
     *
     *  This function requires that the coefficient array
     *  \f$\mathbf{\hat{u}}\f$ provided as \a inarray.
     *
     *  The resulting array
     *  \f$\mathbf{u}[m]=u(\mathbf{\xi}_m)\f$ containing the
     *  expansion evaluated at the quadrature points, is stored
     *  in the \a outarray.
     *
     *  \param inarray contains the values of the expansion
     *  coefficients (input of the function)
     *
     *  \param outarray contains the values of the expansion evaluated
     *  at the quadrature points (output of the function)
     */
    void BwdTrans(const Array<OneD, const NekDouble> &inarray,
                  Array<OneD, NekDouble> &outarray)
    {
        v_BwdTrans(inarray, outarray);
    }
 
    /**
     * @brief This function performs the Forward transformation from
     * physical space to coefficient space.
     */
    inline void FwdTrans(const Array<OneD, const NekDouble> &inarray,
                         Array<OneD, NekDouble> &outarray);
 
    void FwdTransBndConstrained(const Array<OneD, const NekDouble> &inarray,
                                Array<OneD, NekDouble> &outarray)
    {
        v_FwdTransBndConstrained(inarray, outarray);
    }
 
    /** \brief This function integrates the specified function over the
     *  domain
     *
     *  This function is a wrapper around the virtual function
     *  \a v_Integral()
     *
     *  Based on the values of the function evaluated at the quadrature
     *  points (which are stored in \a inarray), this function calculates
     *  the integral of this function over the domain.  This is
     *  equivalent to the numerical evaluation of the operation
     *  \f[ I=\int u(\mathbf{\xi})d \mathbf{\xi}\f]
     *
     *  \param inarray values of the function to be integrated evaluated
     *  at the quadrature points (i.e.
     *  \a inarray[m]=\f$u(\mathbf{\xi}_m)\f$)
     *  \return returns the value of the calculated integral
     *
     *              Inputs:\n
 
    - \a inarray: definition of function to be returned at quadrature point
    of expansion.
 
    Outputs:\n
 
    - returns \f$\int^1_{-1}\int^1_{-1} u(\xi_1, \xi_2) J[i,j] d
    \xi_1 d \xi_2 \f$ where \f$inarray[i,j] =
    u(\xi_{1i},\xi_{2j}) \f$ and \f$ J[i,j] \f$ is the
    Jacobian evaluated at the quadrature point.
     *
     */
    NekDouble Integral(const Array<OneD, const NekDouble> &inarray)
    {
        return v_Integral(inarray);
    }
 
    /** \brief This function fills the array \a outarray with the
     *  \a mode-th mode of the expansion
     *
     *  This function is a wrapper around the virtual function
     *  \a v_FillMode()
     *
     *  The requested mode is evaluated at the quadrature points
     *
     *  \param mode the mode that should be filled
     *  \param outarray contains the values of the \a mode-th mode of the
     *  expansion evaluated at the quadrature points (output of the
     *  function)
     */
    void FillMode(const int mode, Array<OneD, NekDouble> &outarray)
    {
        v_FillMode(mode, outarray);
    }
 
    /** \brief this function calculates the inner product of a given
     *  function \a f with the different modes of the expansion
     *
     *  This function is a wrapper around the virtual function
     *  \a v_IProductWRTBase()
     *
     *  This is equivalent to the numerical evaluation of
     *  \f[ I[p] = \int \phi_p(\mathbf{x}) f(\mathbf{x}) d\mathbf{x}\f]
     *            \f$ \begin{array}{rcl} I_{pq} = (\phi_q \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})
    J_{i,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}
    J_{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$
     *
     *  \param inarray contains the values of the function \a f
     *  evaluated at the quadrature points
     *  \param outarray contains the values of the inner product of \a f
     *  with the different modes, i.e. \f$ outarray[p] = I[p]\f$
     *  (output of the function)
     */
    void IProductWRTBase(const Array<OneD, const NekDouble> &inarray,
                         Array<OneD, NekDouble> &outarray)
    {
        v_IProductWRTBase(inarray, outarray);
    }
 
    void IProductWRTBase(const Array<OneD, const NekDouble> &base,
                         const Array<OneD, const NekDouble> &inarray,
                         Array<OneD, NekDouble> &outarray, int coll_check)
    {
        v_IProductWRTBase(base, inarray, outarray, coll_check);
    }
 
    void IProductWRTDerivBase(const int dir,
                              const Array<OneD, const NekDouble> &inarray,
                              Array<OneD, NekDouble> &outarray)
    {
        v_IProductWRTDerivBase(dir, inarray, outarray);
    }
 
    void IProductWRTDirectionalDerivBase(
        const Array<OneD, const NekDouble> &direction,
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray)
    {
        v_IProductWRTDirectionalDerivBase(direction, inarray, outarray);
    }
 
    /// \brief Get the element id of this expansion when used
    /// in a list by returning value of #m_elmt_id
    inline int GetElmtId()
    {
        return m_elmt_id;
    }
 
    /// \brief Set the element id of this expansion when used
    /// in a list by returning value of #m_elmt_id
    inline void SetElmtId(const int id)
    {
        m_elmt_id = id;
    }
 
    /** \brief this function returns the physical coordinates of the
     *  quadrature points of the expansion
     *
     *  This function is a wrapper around the virtual function
     *  \a v_GetCoords()
     *
     *  \param coords an array containing the coordinates of the
     *  quadrature points (output of the function)
     */
    void GetCoords(Array<OneD, NekDouble> &coords_1,
                   Array<OneD, NekDouble> &coords_2 = NullNekDouble1DArray,
                   Array<OneD, NekDouble> &coords_3 = NullNekDouble1DArray)
    {
        v_GetCoords(coords_1, coords_2, coords_3);
    }
 
    /** \brief given the coordinates of a point of the element in the
     *  local collapsed coordinate system, this function calculates the
     *  physical coordinates of the point
     *
     *  This function is a wrapper around the virtual function
     *  \a v_GetCoord()
     *
     *  \param Lcoords the coordinates in the local collapsed
     *  coordinate system
     *  \param coords the physical coordinates (output of the function)
     */
    void GetCoord(const Array<OneD, const NekDouble> &Lcoord,
                  Array<OneD, NekDouble> &coord)
    {
        v_GetCoord(Lcoord, coord);
    }
 
    inline DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
    {
        return m_stdMatrixManager[mkey];
    }
 
    inline DNekBlkMatSharedPtr GetStdStaticCondMatrix(const StdMatrixKey &mkey)
    {
        return m_stdStaticCondMatrixManager[mkey];
    }
 
    void NormVectorIProductWRTBase(const Array<OneD, const NekDouble> &Fx,
                                   Array<OneD, NekDouble> &outarray)
    {
        v_NormVectorIProductWRTBase(Fx, outarray);
    }
 
    void NormVectorIProductWRTBase(const Array<OneD, const NekDouble> &Fx,
                                   const Array<OneD, NekDouble> &Fy,
                                   Array<OneD, NekDouble> &outarray)
    {
        v_NormVectorIProductWRTBase(Fx, Fy, outarray);
    }
 
    void NormVectorIProductWRTBase(const Array<OneD, const NekDouble> &Fx,
                                   const Array<OneD, const NekDouble> &Fy,
                                   const Array<OneD, const NekDouble> &Fz,
                                   Array<OneD, NekDouble> &outarray)
    {
        v_NormVectorIProductWRTBase(Fx, Fy, Fz, outarray);
    }
 
    void NormVectorIProductWRTBase(
        const Array<OneD, const Array<OneD, NekDouble>> &Fvec,
        Array<OneD, NekDouble> &outarray)
    {
        v_NormVectorIProductWRTBase(Fvec, outarray);
    }
 
    DNekScalBlkMatSharedPtr GetLocStaticCondMatrix(
        const LocalRegions::MatrixKey &mkey)
    {
        return v_GetLocStaticCondMatrix(mkey);
    }
 
    STD_REGIONS_EXPORT void DropLocStaticCondMatrix(
        const LocalRegions::MatrixKey &mkey)
    {
        return v_DropLocStaticCondMatrix(mkey);
    }
 
    int CalcNumberOfCoefficients(const std::vector<unsigned int> &nummodes,
                                 int &modes_offset)
    {
        return v_CalcNumberOfCoefficients(nummodes, modes_offset);
    }
 
    // virtual functions related to LocalRegions
    STD_REGIONS_EXPORT NekDouble
    StdPhysEvaluate(const Array<OneD, const NekDouble> &Lcoord,
                    const Array<OneD, const NekDouble> &physvals);
 
    int GetCoordim()
    {
        return v_GetCoordim();
    }
 
    void GetBoundaryMap(Array<OneD, unsigned int> &outarray)
    {
        v_GetBoundaryMap(outarray);
    }
 
    void GetInteriorMap(Array<OneD, unsigned int> &outarray)
    {
        v_GetInteriorMap(outarray);
    }
 
    int GetVertexMap(const int localVertexId, bool useCoeffPacking = false)
    {
        return v_GetVertexMap(localVertexId, useCoeffPacking);
    }
 
    void GetTraceToElementMap(const int tid,
                              Array<OneD, unsigned int> &maparray,
                              Array<OneD, int> &signarray,
                              Orientation traceOrient = eForwards, int P = -1,
                              int Q = -1)
    {
        v_GetTraceToElementMap(tid, maparray, signarray, traceOrient, P, Q);
    }
 
    void GetTraceCoeffMap(const unsigned int traceid,
                          Array<OneD, unsigned int> &maparray)
    {
        v_GetTraceCoeffMap(traceid, maparray);
    }
 
    void GetElmtTraceToTraceMap(const unsigned int tid,
                                Array<OneD, unsigned int> &maparray,
                                Array<OneD, int> &signarray,
                                Orientation traceOrient = eForwards, int P = -1,
                                int Q = -1)
    {
        v_GetElmtTraceToTraceMap(tid, maparray, signarray, traceOrient, P, Q);
    }
 
    void GetTraceInteriorToElementMap(const int tid,
                                      Array<OneD, unsigned int> &maparray,
                                      Array<OneD, int> &signarray,
                                      const Orientation traceOrient = eForwards)
    {
        v_GetTraceInteriorToElementMap(tid, maparray, signarray, traceOrient);
    }
 
    void GetTraceNumModes(
        const int tid, int &numModes0, int &numModes1,
        const Orientation traceOrient = eDir1FwdDir1_Dir2FwdDir2)
    {
        v_GetTraceNumModes(tid, numModes0, numModes1, traceOrient);
    }
 
    void MultiplyByQuadratureMetric(const Array<OneD, const NekDouble> &inarray,
                                    Array<OneD, NekDouble> &outarray)
    {
        v_MultiplyByQuadratureMetric(inarray, outarray);
    }
 
    void MultiplyByStdQuadratureMetric(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray)
    {
        v_MultiplyByStdQuadratureMetric(inarray, outarray);
    }
 
    // Matrix Routines
 
    /** \brief this function generates the mass matrix
     *  \f$\mathbf{M}[i][j] =
     *  \int \phi_i(\mathbf{x}) \phi_j(\mathbf{x}) d\mathbf{x}\f$
     *
     *  \return returns the mass matrix
     */
 
    STD_REGIONS_EXPORT DNekMatSharedPtr
    CreateGeneralMatrix(const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT void GeneralMatrixOp(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey);
 
    void MassMatrixOp(const Array<OneD, const NekDouble> &inarray,
                      Array<OneD, NekDouble> &outarray,
                      const StdMatrixKey &mkey)
    {
        v_MassMatrixOp(inarray, outarray, mkey);
    }
 
    void LaplacianMatrixOp(const Array<OneD, const NekDouble> &inarray,
                           Array<OneD, NekDouble> &outarray,
                           const StdMatrixKey &mkey)
    {
        v_LaplacianMatrixOp(inarray, outarray, mkey);
    }
 
    void ReduceOrderCoeffs(int numMin,
                           const Array<OneD, const NekDouble> &inarray,
                           Array<OneD, NekDouble> &outarray)
    {
        v_ReduceOrderCoeffs(numMin, inarray, outarray);
    }
 
    void SVVLaplacianFilter(Array<OneD, NekDouble> &array,
                            const StdMatrixKey &mkey)
    {
        v_SVVLaplacianFilter(array, mkey);
    }
 
    void ExponentialFilter(Array<OneD, NekDouble> &array, const NekDouble alpha,
                           const NekDouble exponent, const NekDouble cutoff)
    {
        v_ExponentialFilter(array, alpha, exponent, cutoff);
    }
 
    void LaplacianMatrixOp(const int k1, const int k2,
                           const Array<OneD, const NekDouble> &inarray,
                           Array<OneD, NekDouble> &outarray,
                           const StdMatrixKey &mkey)
    {
        v_LaplacianMatrixOp(k1, k2, inarray, outarray, mkey);
    }
 
    void WeakDerivMatrixOp(const int i,
                           const Array<OneD, const NekDouble> &inarray,
                           Array<OneD, NekDouble> &outarray,
                           const StdMatrixKey &mkey)
    {
        v_WeakDerivMatrixOp(i, inarray, outarray, mkey);
    }
 
    void WeakDirectionalDerivMatrixOp(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
    {
        v_WeakDirectionalDerivMatrixOp(inarray, outarray, mkey);
    }
 
    void MassLevelCurvatureMatrixOp(const Array<OneD, const NekDouble> &inarray,
                                    Array<OneD, NekDouble> &outarray,
                                    const StdMatrixKey &mkey)
    {
        v_MassLevelCurvatureMatrixOp(inarray, outarray, mkey);
    }
 
    void LinearAdvectionDiffusionReactionMatrixOp(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey,
        bool addDiffusionTerm = true)
    {
        v_LinearAdvectionDiffusionReactionMatrixOp(inarray, outarray, mkey,
                                                   addDiffusionTerm);
    }
 
    /**
     * @param   inarray     Input array @f$ \mathbf{u} @f$.
     * @param   outarray    Output array @f$ \boldsymbol{\nabla^2u}
     *                          + \lambda \boldsymbol{u} @f$.
     * @param   mkey
     */
    void HelmholtzMatrixOp(const Array<OneD, const NekDouble> &inarray,
                           Array<OneD, NekDouble> &outarray,
                           const StdMatrixKey &mkey)
    {
        v_HelmholtzMatrixOp(inarray, outarray, mkey);
    }
 
    DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
    {
        return v_GenMatrix(mkey);
    }
 
    void PhysDeriv(const Array<OneD, const NekDouble> &inarray,
                   Array<OneD, NekDouble> &out_d0,
                   Array<OneD, NekDouble> &out_d1 = NullNekDouble1DArray,
                   Array<OneD, NekDouble> &out_d2 = NullNekDouble1DArray)
    {
        v_PhysDeriv(inarray, out_d0, out_d1, out_d2);
    }
 
    void PhysDeriv(const int dir, const Array<OneD, const NekDouble> &inarray,
                   Array<OneD, NekDouble> &outarray)
    {
        v_PhysDeriv(dir, inarray, outarray);
    }
 
    void PhysDeriv_s(const Array<OneD, const NekDouble> &inarray,
                     Array<OneD, NekDouble> &out_ds)
    {
        v_PhysDeriv_s(inarray, out_ds);
    }
 
    void PhysDeriv_n(const Array<OneD, const NekDouble> &inarray,
                     Array<OneD, NekDouble> &out_dn)
    {
        v_PhysDeriv_n(inarray, out_dn);
    }
 
    void PhysDirectionalDeriv(const Array<OneD, const NekDouble> &inarray,
                              const Array<OneD, const NekDouble> &direction,
                              Array<OneD, NekDouble> &outarray)
    {
        v_PhysDirectionalDeriv(inarray, direction, outarray);
    }
 
    void StdPhysDeriv(const Array<OneD, const NekDouble> &inarray,
                      Array<OneD, NekDouble> &out_d0,
                      Array<OneD, NekDouble> &out_d1 = NullNekDouble1DArray,
                      Array<OneD, NekDouble> &out_d2 = NullNekDouble1DArray)
    {
        v_StdPhysDeriv(inarray, out_d0, out_d1, out_d2);
    }
 
    void StdPhysDeriv(const int dir,
                      const Array<OneD, const NekDouble> &inarray,
                      Array<OneD, NekDouble> &outarray)
    {
        v_StdPhysDeriv(dir, inarray, outarray);
    }
 
    /** \brief This function evaluates the expansion at a single
     *  (arbitrary) point of the domain
     *
     *  This function is a wrapper around the virtual function
     *  \a v_PhysEvaluate()
     *
     *  Based on the value of the expansion at the quadrature
     *  points provided in \a physvals, this function
     *  calculates the value of the expansion at an arbitrary
     *  single points (with coordinates \f$ \mathbf{x_c}\f$
     *  given by the pointer \a coords). This operation,
     *  equivalent to \f[ u(\mathbf{x_c}) = \sum_p
     *  \phi_p(\mathbf{x_c}) \hat{u}_p \f] is evaluated using
     *  Lagrangian interpolants through the quadrature points:
     *  \f[ u(\mathbf{x_c}) = \sum_p h_p(\mathbf{x_c}) u_p\f]
     *
     *  \param coords the coordinates of the single point
     *  \param physvals the interpolated field at the quadrature points
     *
     *  \return returns the value of the expansion at the
     *  single point
     */
    NekDouble PhysEvaluate(const Array<OneD, const NekDouble> &coords,
                           const Array<OneD, const NekDouble> &physvals)
    {
        return v_PhysEvaluate(coords, physvals);
    }
 
    /** \brief This function evaluates the first derivative of the expansion
     * at a single (arbitrary) point of the domain
     *
     *  This function is a wrapper around the virtual function
     *  \a v_PhysEvaluate()
 
     *  Based on the value of the expansion at the quadrature
     *  points provided in \a physvals, this function
     *  calculates the value of the expansion at a set of points
     * given in \a coords
     */
    inline NekDouble PhysEvaluate(const Array<OneD, NekDouble> &coord,
                                  const Array<OneD, const NekDouble> &inarray,
                                  std::array<NekDouble, 3> &firstOrderDerivs)
 
    {
        return v_PhysEvaluate(coord, inarray, firstOrderDerivs);
    }
 
    inline NekDouble PhysEvaluate(const Array<OneD, NekDouble> &coord,
                                  const Array<OneD, const NekDouble> &inarray,
                                  std::array<NekDouble, 3> &firstOrderDerivs,
                                  std::array<NekDouble, 6> &secondOrderDerivs)
 
    {
        return v_PhysEvaluate(coord, inarray, firstOrderDerivs,
                              secondOrderDerivs);
    }
 
    /** \brief This function evaluates the expansion at a single
     *  (arbitrary) point of the domain
     *
     *  This function is a wrapper around the virtual function
     *  \a v_PhysEvaluate()
     *
     *  Based on the value of the expansion at the quadrature
     *  points provided in \a physvals, this function
     *  calculates the value of the expansion at an arbitrary
     *  single points associated with the interpolation
     *  matrices provided in \f$ I \f$.
     *
     *  \param I an Array of lagrange interpolantes evaluated
     *  at the coordinate and going through the local physical
     *  quadrature
     *  \param physvals the interpolated field at the quadrature points
     *
     *  \return returns the value of the expansion at the
     *  single point
     */
    NekDouble PhysEvaluate(const Array<OneD, DNekMatSharedPtr> &I,
                           const Array<OneD, const NekDouble> &physvals)
    {
        return v_PhysEvaluate(I, physvals);
    }
 
    /**
     * @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 PhysEvaluateBasis(const Array<OneD, const NekDouble> &coords,
                                int mode)
    {
        return v_PhysEvaluateBasis(coords, mode);
    }
 
    /**
     * \brief Convert local cartesian coordinate \a xi into local
     * collapsed coordinates \a eta
     **/
    void LocCoordToLocCollapsed(const Array<OneD, const NekDouble> &xi,
                                Array<OneD, NekDouble> &eta)
    {
        v_LocCoordToLocCollapsed(xi, eta);
    }
 
    /**
     * \brief Convert local collapsed coordinates \a eta into local
     * cartesian coordinate \a xi
     **/
    void LocCollapsedToLocCoord(const Array<OneD, const NekDouble> &eta,
                                Array<OneD, NekDouble> &xi)
    {
        v_LocCollapsedToLocCoord(eta, xi);
    }
 
    STD_REGIONS_EXPORT virtual int v_CalcNumberOfCoefficients(
        const std::vector<unsigned int> &nummodes, int &modes_offset);
 
    STD_REGIONS_EXPORT virtual void v_NormVectorIProductWRTBase(
        const Array<OneD, const NekDouble> &Fx,
        Array<OneD, NekDouble> &outarray);
 
    STD_REGIONS_EXPORT virtual void v_NormVectorIProductWRTBase(
        const Array<OneD, const NekDouble> &Fx,
        const Array<OneD, const NekDouble> &Fy,
        Array<OneD, NekDouble> &outarray);
 
    STD_REGIONS_EXPORT virtual void v_NormVectorIProductWRTBase(
        const Array<OneD, const NekDouble> &Fx,
        const Array<OneD, const NekDouble> &Fy,
        const Array<OneD, const NekDouble> &Fz,
        Array<OneD, NekDouble> &outarray);
 
    STD_REGIONS_EXPORT virtual void v_NormVectorIProductWRTBase(
        const Array<OneD, const Array<OneD, NekDouble>> &Fvec,
        Array<OneD, NekDouble> &outarray);
 
    STD_REGIONS_EXPORT virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(
        const LocalRegions::MatrixKey &mkey);
 
    STD_REGIONS_EXPORT virtual void v_DropLocStaticCondMatrix(
        const LocalRegions::MatrixKey &mkey);
 
    /** \brief Function to evaluate the discrete \f$ L_\infty\f$
     *  error \f$ |\epsilon|_\infty = \max |u - u_{exact}|\f$ where \f$
     *    u_{exact}\f$ is given by the array \a sol.
     *
     *    This function takes the physical value space array \a m_phys as
     *  approximate solution
     *
     *  \param sol array of solution function  at physical quadrature
     *  points
     *  \return returns the \f$ L_\infty \f$ error as a NekDouble.
     */
    STD_REGIONS_EXPORT NekDouble
    Linf(const Array<OneD, const NekDouble> &phys,
         const Array<OneD, const NekDouble> &sol = NullNekDouble1DArray);
 
    /** \brief Function to evaluate the discrete \f$ L_2\f$ error,
     *  \f$ | \epsilon |_{2} = \left [ \int^1_{-1} [u - u_{exact}]^2
     *  dx \right]^{1/2} d\xi_1 \f$ where \f$ u_{exact}\f$ is given by
     *  the array \a sol.
     *
     *    This function takes the physical value space array \a m_phys as
     *  approximate solution
     *
     *  \param sol array of solution function  at physical quadrature
     *  points
     *  \return returns the \f$ L_2 \f$ error as a double.
     */
    STD_REGIONS_EXPORT NekDouble
    L2(const Array<OneD, const NekDouble> &phys,
       const Array<OneD, const NekDouble> &sol = NullNekDouble1DArray);
 
    /** \brief Function to evaluate the discrete \f$ H^1\f$
     *  error, \f$ | \epsilon |^1_{2} = \left [ \int^1_{-1} [u -
     *  u_{exact}]^2 + \nabla(u - u_{exact})\cdot\nabla(u -
     *  u_{exact})\cdot dx \right]^{1/2} d\xi_1 \f$ where \f$
     *  u_{exact}\f$ is given by the array \a sol.
     *
     *    This function takes the physical value space array
     *  \a m_phys as approximate solution
     *
     *  \param sol array of solution function  at physical quadrature
     *  points
     *  \return returns the \f$ H_1 \f$ error as a double.
     */
    STD_REGIONS_EXPORT NekDouble
    H1(const Array<OneD, const NekDouble> &phys,
       const Array<OneD, const NekDouble> &sol = NullNekDouble1DArray);
 
    // I/O routines
    const LibUtilities::PointsKeyVector GetPointsKeys() const
    {
        LibUtilities::PointsKeyVector p;
        p.reserve(m_base.size());
        for (size_t i = 0; i < m_base.size(); ++i)
        {
            p.push_back(m_base[i]->GetPointsKey());
        }
        return p;
    }
 
    STD_REGIONS_EXPORT DNekMatSharedPtr BuildInverseTransformationMatrix(
        const DNekScalMatSharedPtr &m_transformationmatrix)
    {
        return v_BuildInverseTransformationMatrix(m_transformationmatrix);
    }
 
    /** \brief This function performs an interpolation from
     * the physical space points provided at input into an
     * array of equispaced points which are not the collapsed
     * coordinate. So for a tetrahedron you will only get a
     * tetrahedral number of values.
     *
     * This is primarily used for output purposes to get a
     * better distribution of points more suitable for most
     * postprocessing
     */
    STD_REGIONS_EXPORT void PhysInterpToSimplexEquiSpaced(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, int npset = -1);
 
    /** \brief This function provides the connectivity of
     *   local simplices (triangles or tets) to connect the
     *   equispaced data points provided by
     *   PhysInterpToSimplexEquiSpaced
     *
     *  This is a virtual call to the function
     *  \a v_GetSimplexEquiSpaceConnectivity
     */
    STD_REGIONS_EXPORT void GetSimplexEquiSpacedConnectivity(
        Array<OneD, int> &conn, bool standard = true)
    {
        v_GetSimplexEquiSpacedConnectivity(conn, standard);
    }
 
    /** \brief This function performs a
     * projection/interpolation from the equispaced points
     * sometimes used in post-processing onto the coefficient
     * space
     *
     * This is primarily used for output purposes to use a
     * more even distribution of points more suitable for alot of
     * postprocessing
     */
    STD_REGIONS_EXPORT void EquiSpacedToCoeffs(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray);
 
    template <class T> std::shared_ptr<T> as()
    {
        return std::dynamic_pointer_cast<T>(shared_from_this());
    }
 
    void IProductWRTBase_SumFac(const Array<OneD, const NekDouble> &inarray,
                                Array<OneD, NekDouble> &outarray,
                                bool multiplybyweights = true)
    {
        v_IProductWRTBase_SumFac(inarray, outarray, multiplybyweights);
    }
 
    STD_REGIONS_EXPORT void GenStdMatBwdDeriv(const int dir,
                                              DNekMatSharedPtr &mat)
    {
        v_GenStdMatBwdDeriv(dir, mat);
    }
 
protected:
    Array<OneD, LibUtilities::BasisSharedPtr>
        m_base; /**< Bases needed for the expansion */
    int m_elmt_id;
    int m_ncoeffs; /**< Total number of coefficients used in the expansion */
 
    LibUtilities::NekManager<StdMatrixKey, DNekMat, StdMatrixKey::opLess>
        m_stdMatrixManager;
    LibUtilities::NekManager<StdMatrixKey, DNekBlkMat, StdMatrixKey::opLess>
        m_stdStaticCondMatrixManager;
 
    DNekMatSharedPtr CreateStdMatrix(const StdMatrixKey &mkey)
    {
        return v_CreateStdMatrix(mkey);
    }
 
    /** \brief Create the static condensation of a matrix when
        using a boundary interior decomposition
 
        If a matrix system can be represented by
        \f$ Mat = \left [ \begin{array}{cc}
        A & B \\
        C & D \end{array} \right ] \f$
        This routine creates a matrix containing the statically
        condense system of the form
        \f$ Mat = \left [ \begin{array}{cc}
        A - B D^{-1} C & B D^{-1} \\
        D^{-1} C       & D^{-1} \end{array} \right ] \f$
    **/
    STD_REGIONS_EXPORT DNekBlkMatSharedPtr
    CreateStdStaticCondMatrix(const StdMatrixKey &mkey);
 
    void BwdTrans_SumFac(const Array<OneD, const NekDouble> &inarray,
                         Array<OneD, NekDouble> &outarray)
    {
        v_BwdTrans_SumFac(inarray, outarray);
    }
 
    void IProductWRTDerivBase_SumFac(
        const int dir, const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray)
    {
        v_IProductWRTDerivBase_SumFac(dir, inarray, outarray);
    }
 
    void IProductWRTDirectionalDerivBase_SumFac(
        const Array<OneD, const NekDouble> &direction,
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray)
    {
        v_IProductWRTDirectionalDerivBase_SumFac(direction, inarray, outarray);
    }
 
    // The term _MatFree denotes that the action of the
    // MatrixOperation is done withouth actually using the
    // matrix (which then needs to be stored/calculated).
    // Although this does not strictly mean that no matrix
    // operations are involved in the evaluation of the
    // operation, we use this term in the same context used as
    // in the following paper: R. C. Kirby, M. G. Knepley,
    // A. Logg, and L. R. Scott, "Optimizing the evaluation of
    // finite element matrices," SISC 27:741-758 (2005)
    STD_REGIONS_EXPORT void GeneralMatrixOp_MatFree(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT void MassMatrixOp_MatFree(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey);
 
    void LaplacianMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
                                   Array<OneD, NekDouble> &outarray,
                                   const StdMatrixKey &mkey)
    {
        v_LaplacianMatrixOp_MatFree(inarray, outarray, mkey);
    }
 
    STD_REGIONS_EXPORT void LaplacianMatrixOp_MatFree_Kernel(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, Array<OneD, NekDouble> &wsp)
    {
        v_LaplacianMatrixOp_MatFree_Kernel(inarray, outarray, wsp);
    }
 
    STD_REGIONS_EXPORT void LaplacianMatrixOp_MatFree_GenericImpl(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT void LaplacianMatrixOp_MatFree(
        const int k1, const int k2, const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT void WeakDerivMatrixOp_MatFree(
        const int i, const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT void WeakDirectionalDerivMatrixOp_MatFree(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT void MassLevelCurvatureMatrixOp_MatFree(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT void LinearAdvectionDiffusionReactionMatrixOp_MatFree(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey,
        bool addDiffusionTerm = true);
 
    void HelmholtzMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
                                   Array<OneD, NekDouble> &outarray,
                                   const StdMatrixKey &mkey)
    {
        v_HelmholtzMatrixOp_MatFree(inarray, outarray, mkey);
    }
 
    STD_REGIONS_EXPORT void HelmholtzMatrixOp_MatFree_GenericImpl(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT virtual void v_SetCoeffsToOrientation(
        StdRegions::Orientation dir, Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray);
 
    STD_REGIONS_EXPORT virtual void v_SetCoeffsToOrientation(
        Array<OneD, NekDouble> &coeffs, StdRegions::Orientation dir);
 
    STD_REGIONS_EXPORT virtual NekDouble v_StdPhysEvaluate(
        const Array<OneD, const NekDouble> &Lcoord,
        const Array<OneD, const NekDouble> &physvals);
 
    STD_REGIONS_EXPORT virtual void v_GenStdMatBwdDeriv(const int dir,
                                                        DNekMatSharedPtr &mat)
    {
        boost::ignore_unused(dir, mat);
        NEKERROR(ErrorUtil::efatal, "not defined");
    }
 
    STD_REGIONS_EXPORT virtual void v_MultiplyByStdQuadratureMetric(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray);
 
    /**
     * @brief This function performs the barycentric interpolation of
     * the polynomial stored in @p coord at a point @p physvals using
     * barycentric interpolation weights in direction @tparam DIR.
     * It can also perform the barycentric interpolation of the
     * derivative of the polynomial if @tparam DERIV is set to true,
     * which outputs in to @param deriv.
     *
     * This method is intended to be used a helper function for
     * StdExpansion::PhysEvaluate and its elemental instances, so that
     * the calling method should provide @p coord for x, y and z
     * sequentially and the appropriate @p physvals and @p weights for
     * that particular direction.
     *
     * @param  coord    The coordinate of the single point.
     * @param  physvals The polynomial stored at each quadrature point.
     * @param  deriv    The value of the derivative.
     * @param  deriv    The value of the 2nd derivative.
     * @tparam DIR      The direction of evaluation.
     * @tparam DERIV    Bool to find derivative.
     *
     * @return The value of @p physvals at @p coord in direction @p dir.
     */
    template <int DIR, bool DERIV = false, bool DERIV2 = false>
    inline NekDouble BaryEvaluate(const NekDouble &coord,
                                  const NekDouble *physvals, NekDouble &deriv,
                                  NekDouble &deriv2)
    {
        NekDouble numer1 = 0.0, numer2 = 0.0, numer3 = 0.0, numer4 = 0.0,
                  numer5 = 0.0, denom = 0.0;
 
        ASSERTL2(DIR < m_base.size(),
                 "Direction should be less than shape dimension.");
 
        const Array<OneD, const NekDouble> &z  = m_base[DIR]->GetZ();
        const Array<OneD, const NekDouble> &bw = m_base[DIR]->GetBaryWeights();
 
        const size_t nquad = z.size();
 
        for (size_t i = 0; i < nquad; ++i)
        {
            NekDouble xdiff = z[i] - coord;
            NekDouble pval  = physvals[i];
 
            /*
             * (in this specific case) you actually
             * want to do the comparison exactly
             * (believe it or not!) See chapter 7 of
             * the paper here:
             *https://people.maths.ox.ac.uk/trefethen/barycentric.pdf
             */
            if ((!DERIV && xdiff == 0.0) ||
                ((DERIV || DERIV2) && std::abs(xdiff) < 1e-15))
            {
 
                if (DERIV2)
                {
                    DNekMatSharedPtr D0 = m_base[DIR]->GetD();
 
                    // take ith row of z and multiply with physvals
                    Array<OneD, NekDouble> tmp(nquad);
                    for (int kk = 0; kk < nquad; kk++)
                    {
                        tmp[kk] = Vmath::Dot(nquad, &(D0->GetPtr())[kk], nquad,
                                             &physvals[0], 1);
                    }
 
                    deriv2 = Vmath::Dot(nquad, &(D0->GetPtr())[i], nquad,
                                        &tmp[0], 1);
                    deriv  = tmp[i];
                }
                else if (DERIV)
                {
                    DNekMatSharedPtr D0 = m_base[DIR]->GetD();
 
                    // take ith row of z and multiply with physvals
                    deriv = Vmath::Dot(z.size(), &(D0->GetPtr())[i], z.size(),
                                       &physvals[0], 1);
                }
 
                return pval;
            }
 
            NekDouble tmp = bw[i] / xdiff;
            numer1 += tmp * pval;
            denom += tmp;
 
            if (DERIV || DERIV2)
            {
                NekDouble tmp2 = tmp / xdiff;
                numer2 += tmp2 * pval;
                numer3 += tmp2;
 
                if (DERIV2)
                {
                    NekDouble tmp3 = tmp2 / xdiff;
                    numer4 += tmp3 * pval;
                    numer5 += tmp3;
                }
            }
        }
 
        if (DERIV || DERIV2)
        {
            NekDouble denomdenom   = denom * denom;
            NekDouble numer1numer3 = numer1 * numer3;
 
            deriv = (numer2 * denom - numer1numer3) / (denomdenom);
 
            if (DERIV2)
            {
                deriv2 = (2.0 * numer4 / denom) -
                         (2.0 * numer5 * numer1) / (denomdenom) -
                         (2.0 * numer2 * numer3) / (denomdenom) +
                         (2.0 * numer3 * numer1numer3) / (denomdenom * denom);
            }
        }
 
        return numer1 / denom;
    }
 
    /**
     * @brief This function evaluates the basis function mode @p mode at
     * a point @p coords of the domain in direction @dir.
     *
     * @param coord   The coordinate inside the standard region.
     * @param mode    The mode number to be evaluated of #m_base[dir]
     * @param dir     The direction of interpolation.
     *
     * @return The value of the basis function @p mode at @p coords in
     *         direction @p dir.
     */
    template <int DIR>
    inline NekDouble BaryEvaluateBasis(const NekDouble &coord, const int &mode)
    {
        const int nquad = m_base[DIR]->GetNumPoints();
        return BaryEvaluate<DIR>(coord,
                                 &(m_base[DIR]->GetBdata())[0] + nquad * mode);
    }
 
    /**
     * @brief Helper function to pass an unused value by reference into
     * BaryEvaluate.
     *
     * @param  coord    The coordinate of the single point.
     * @param  physvals The polynomial stored at each quadrature point.
     * @tparam DIR      The direction of evaluation.
     * @tparam DERIV    Bool to find derivative.
     *
     * @return The value of @p physvals at @p coord in direction @p dir.
     */
    template <int DIR, bool DERIV = false, bool DERIV2 = false>
    inline NekDouble BaryEvaluate(const NekDouble &coord,
                                  const NekDouble *physvals)
    {
        NekDouble unusedValue = 0.0;
        return BaryEvaluate<DIR, DERIV, DERIV2>(coord, physvals, unusedValue,
                                                unusedValue);
    }
 
    template <int DIR, bool DERIV = false, bool DERIV2 = false>
    inline NekDouble BaryEvaluate(const NekDouble &coord,
                                  const NekDouble *physvals, NekDouble &deriv)
    {
        NekDouble unusedValue = 0.0;
        return BaryEvaluate<DIR, DERIV, DERIV2>(coord, physvals, deriv,
                                                unusedValue);
    }
 
private:
    // Virtual functions
    STD_REGIONS_EXPORT virtual int v_GetNverts() const  = 0;
    STD_REGIONS_EXPORT virtual int v_GetNtraces() const = 0;
 
    STD_REGIONS_EXPORT virtual int v_NumBndryCoeffs() const   = 0;
    STD_REGIONS_EXPORT virtual int v_NumDGBndryCoeffs() const = 0;
 
    STD_REGIONS_EXPORT virtual int v_GetTraceNcoeffs(const int i) const    = 0;
    STD_REGIONS_EXPORT virtual int v_GetTraceIntNcoeffs(const int i) const = 0;
    STD_REGIONS_EXPORT virtual int v_GetTraceNumPoints(const int i) const  = 0;
 
    STD_REGIONS_EXPORT virtual const LibUtilities::BasisKey v_GetTraceBasisKey(
        const int i, const int k) const;
 
    STD_REGIONS_EXPORT virtual LibUtilities::PointsKey v_GetTracePointsKey(
        const int i, const int j) const;
 
    STD_REGIONS_EXPORT virtual const LibUtilities::PointsKey v_GetNodalPointsKey()
        const;
 
    STD_REGIONS_EXPORT virtual LibUtilities::ShapeType v_DetShapeType()
        const = 0;
 
    STD_REGIONS_EXPORT virtual std::shared_ptr<StdExpansion> v_GetStdExp()
        const;
 
    STD_REGIONS_EXPORT virtual std::shared_ptr<StdExpansion> v_GetLinStdExp(
        void) const;
 
    STD_REGIONS_EXPORT virtual int v_GetShapeDimension() const = 0;
 
    STD_REGIONS_EXPORT virtual bool v_IsBoundaryInteriorExpansion() const;
 
    STD_REGIONS_EXPORT virtual bool v_IsNodalNonTensorialExp();
 
    STD_REGIONS_EXPORT virtual void v_BwdTrans(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray) = 0;
 
    /**
     * @brief Transform a given function from physical quadrature space
     * to coefficient space.
     * @see StdExpansion::FwdTrans
     */
    STD_REGIONS_EXPORT virtual void v_FwdTrans(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray) = 0;
 
    /**
     * @brief Calculates the inner product of a given function \a f
     * with the different modes of the expansion
     */
    STD_REGIONS_EXPORT virtual void v_IProductWRTBase(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray) = 0;
 
    STD_REGIONS_EXPORT virtual void v_IProductWRTBase(
        const Array<OneD, const NekDouble> &base,
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, int coll_check)
    {
        boost::ignore_unused(base, inarray, outarray, coll_check);
        NEKERROR(ErrorUtil::efatal, "StdExpansion::v_IProductWRTBase has no "
                                    "(and should have no) implementation");
    }
 
    STD_REGIONS_EXPORT virtual void v_IProductWRTDerivBase(
        const int dir, const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray);
 
    STD_REGIONS_EXPORT virtual void v_IProductWRTDirectionalDerivBase(
        const Array<OneD, const NekDouble> &direction,
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray);
 
    STD_REGIONS_EXPORT virtual void v_FwdTransBndConstrained(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray);
 
    STD_REGIONS_EXPORT virtual NekDouble v_Integral(
        const Array<OneD, const NekDouble> &inarray);
 
    STD_REGIONS_EXPORT virtual void v_PhysDeriv(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &out_d1, Array<OneD, NekDouble> &out_d2,
        Array<OneD, NekDouble> &out_d3);
 
    STD_REGIONS_EXPORT virtual void v_PhysDeriv_s(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &out_ds);
 
    STD_REGIONS_EXPORT virtual void v_PhysDeriv_n(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &out_dn);
 
    STD_REGIONS_EXPORT virtual void v_PhysDeriv(
        const int dir, const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &out_d0);
 
    STD_REGIONS_EXPORT virtual void v_PhysDirectionalDeriv(
        const Array<OneD, const NekDouble> &inarray,
        const Array<OneD, const NekDouble> &direction,
        Array<OneD, NekDouble> &outarray);
 
    STD_REGIONS_EXPORT virtual void v_StdPhysDeriv(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &out_d1, Array<OneD, NekDouble> &out_d2,
        Array<OneD, NekDouble> &out_d3);
 
    STD_REGIONS_EXPORT virtual void v_StdPhysDeriv(
        const int dir, const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray);
 
    STD_REGIONS_EXPORT virtual NekDouble v_PhysEvaluate(
        const Array<OneD, const NekDouble> &coords,
        const Array<OneD, const NekDouble> &physvals);
 
    STD_REGIONS_EXPORT virtual NekDouble v_PhysEvaluate(
        const Array<OneD, DNekMatSharedPtr> &I,
        const Array<OneD, const NekDouble> &physvals);
 
    STD_REGIONS_EXPORT virtual NekDouble v_PhysEvaluate(
        const Array<OneD, NekDouble> &coord,
        const Array<OneD, const NekDouble> &inarray,
        std::array<NekDouble, 3> &firstOrderDerivs);
 
    STD_REGIONS_EXPORT virtual NekDouble v_PhysEvaluate(
        const Array<OneD, NekDouble> &coord,
        const Array<OneD, const NekDouble> &inarray,
        std::array<NekDouble, 3> &firstOrderDerivs,
        std::array<NekDouble, 6> &secondOrderDerivs);
 
    STD_REGIONS_EXPORT virtual NekDouble v_PhysEvaluateBasis(
        const Array<OneD, const NekDouble> &coords, int mode);
 
    STD_REGIONS_EXPORT virtual void v_LocCoordToLocCollapsed(
        const Array<OneD, const NekDouble> &xi, Array<OneD, NekDouble> &eta);
 
    STD_REGIONS_EXPORT virtual void v_LocCollapsedToLocCoord(
        const Array<OneD, const NekDouble> &eta, Array<OneD, NekDouble> &xi);
 
    STD_REGIONS_EXPORT virtual void v_FillMode(
        const int mode, Array<OneD, NekDouble> &outarray);
 
    STD_REGIONS_EXPORT virtual DNekMatSharedPtr v_GenMatrix(
        const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT virtual DNekMatSharedPtr v_CreateStdMatrix(
        const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT virtual void v_GetCoords(
        Array<OneD, NekDouble> &coords_0, Array<OneD, NekDouble> &coords_1,
        Array<OneD, NekDouble> &coords_2);
 
    STD_REGIONS_EXPORT virtual void v_GetCoord(
        const Array<OneD, const NekDouble> &Lcoord,
        Array<OneD, NekDouble> &coord);
 
    STD_REGIONS_EXPORT virtual int v_GetCoordim() const;
 
    STD_REGIONS_EXPORT virtual void v_GetBoundaryMap(
        Array<OneD, unsigned int> &outarray);
 
    STD_REGIONS_EXPORT virtual void v_GetInteriorMap(
        Array<OneD, unsigned int> &outarray);
 
    STD_REGIONS_EXPORT virtual int v_GetVertexMap(int localVertexId,
                                                  bool useCoeffPacking = false);
 
    STD_REGIONS_EXPORT virtual void v_GetTraceToElementMap(
        const int tid, Array<OneD, unsigned int> &maparray,
        Array<OneD, int> &signarray, Orientation traceOrient = eForwards,
        int P = -1, int Q = -1);
 
    STD_REGIONS_EXPORT virtual void v_GetTraceCoeffMap(
        const unsigned int traceid, Array<OneD, unsigned int> &maparray);
 
    STD_REGIONS_EXPORT virtual void v_GetElmtTraceToTraceMap(
        const unsigned int tid, Array<OneD, unsigned int> &maparray,
        Array<OneD, int> &signarray, Orientation traceOrient = eForwards,
        int P = -1, int Q = -1);
 
    STD_REGIONS_EXPORT virtual void v_GetTraceInteriorToElementMap(
        const int eid, Array<OneD, unsigned int> &maparray,
        Array<OneD, int> &signarray, const Orientation traceOrient = eForwards);
 
    STD_REGIONS_EXPORT virtual void v_GetTraceNumModes(
        const int fid, int &numModes0, int &numModes1,
        Orientation traceOrient = eDir1FwdDir1_Dir2FwdDir2);
 
    STD_REGIONS_EXPORT virtual void v_GetVertexPhysVals(
        const int vertex, const Array<OneD, const NekDouble> &inarray,
        NekDouble &outarray);
 
    STD_REGIONS_EXPORT virtual void v_MultiplyByQuadratureMetric(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray);
 
    STD_REGIONS_EXPORT virtual void v_BwdTrans_SumFac(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray);
 
    STD_REGIONS_EXPORT virtual void v_IProductWRTBase_SumFac(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, bool multiplybyweights = true);
 
    STD_REGIONS_EXPORT virtual void v_IProductWRTDerivBase_SumFac(
        const int dir, const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray);
 
    STD_REGIONS_EXPORT virtual void v_IProductWRTDirectionalDerivBase_SumFac(
        const Array<OneD, const NekDouble> &direction,
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray);
 
    STD_REGIONS_EXPORT virtual void v_MassMatrixOp(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT virtual void v_LaplacianMatrixOp(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT virtual void v_SVVLaplacianFilter(
        Array<OneD, NekDouble> &array, const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT virtual void v_ExponentialFilter(
        Array<OneD, NekDouble> &array, const NekDouble alpha,
        const NekDouble exponent, const NekDouble cutoff);
 
    STD_REGIONS_EXPORT virtual void v_ReduceOrderCoeffs(
        int numMin, const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray);
 
    STD_REGIONS_EXPORT virtual void v_LaplacianMatrixOp(
        const int k1, const int k2, const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT virtual void v_WeakDerivMatrixOp(
        const int i, const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT virtual void v_WeakDirectionalDerivMatrixOp(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT virtual void v_MassLevelCurvatureMatrixOp(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT virtual void v_LinearAdvectionDiffusionReactionMatrixOp(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey,
        bool addDiffusionTerm = true);
 
    STD_REGIONS_EXPORT virtual void v_HelmholtzMatrixOp(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT virtual void v_LaplacianMatrixOp_MatFree(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT virtual void v_LaplacianMatrixOp_MatFree_Kernel(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, Array<OneD, NekDouble> &wsp);
 
    STD_REGIONS_EXPORT virtual void v_HelmholtzMatrixOp_MatFree(
        const Array<OneD, const NekDouble> &inarray,
        Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey);
 
    STD_REGIONS_EXPORT virtual DNekMatSharedPtr v_BuildInverseTransformationMatrix(
        const DNekScalMatSharedPtr &m_transformationmatrix);
 
    STD_REGIONS_EXPORT virtual void v_GetSimplexEquiSpacedConnectivity(
        Array<OneD, int> &conn, bool standard = true);
};
 
typedef std::shared_ptr<StdExpansion> StdExpansionSharedPtr;
typedef std::vector<StdExpansionSharedPtr> StdExpansionVector;
 
/**
 *  This function is a wrapper around the virtual function
 *  \a v_FwdTrans()
 *
 *  Given a function evaluated at the quadrature points, this
 *  function calculates the expansion coefficients such that the
 *  resulting expansion approximates the original function.
 *
 *  The calculation of the expansion coefficients is done using a
 *  Galerkin projection. This is equivalent to the operation:
 *  \f[ \mathbf{\hat{u}} = \mathbf{M}^{-1} \mathbf{I}\f]
 *  where
 *  - \f$\mathbf{M}[p][q]= \int\phi_p(\mathbf{\xi})\phi_q(
 *  \mathbf{\xi}) d\mathbf{\xi}\f$ is the Mass matrix
 *  - \f$\mathbf{I}[p] = \int\phi_p(\mathbf{\xi}) u(\mathbf{\xi})
 *  d\mathbf{\xi}\f$
 *
 *  This function takes the array \a inarray as the values of the
 *  function evaluated at the quadrature points
 *  (i.e. \f$\mathbf{u}\f$),
 *  and stores the resulting coefficients \f$\mathbf{\hat{u}}\f$
 *  in the \a outarray
 *
 *  @param inarray array of the function discretely evaluated at the
 *  quadrature points
 *
 *  @param outarray array of the function coefficieints
 */
inline void StdExpansion::FwdTrans(const Array<OneD, const NekDouble> &inarray,
                                   Array<OneD, NekDouble> &outarray)
{
    v_FwdTrans(inarray, outarray);
}
 
} // namespace StdRegions
} // namespace Nektar
 
#endif // STANDARDDEXPANSION_H