ExpList.h

///////////////////////////////////////////////////////////////////////////////
//
// File ExpList.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).
//
// License for the specific language governing rights and limitations under
// 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: Expansion list top class definition
//
///////////////////////////////////////////////////////////////////////////////

#ifndef NEKTAR_LIBS_MULTIREGIONS_EXPLIST_H
#define NEKTAR_LIBS_MULTIREGIONS_EXPLIST_H

#include <LibUtilities/Communication/Transposition.h>
#include <LibUtilities/Communication/Comm.h>
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <SpatialDomains/MeshGraph.h>
#include <LocalRegions/Expansion.h>
#include <Collections/Collection.h>
#include <MultiRegions/MultiRegionsDeclspec.h>
#include <MultiRegions/MultiRegions.hpp>
#include <MultiRegions/GlobalMatrix.h>
#include <MultiRegions/GlobalMatrixKey.h>
#include <MultiRegions/GlobalOptimizationParameters.h>
#include <MultiRegions/AssemblyMap/AssemblyMap.h>
#include <boost/enable_shared_from_this.hpp>
#include <tinyxml.h>

namespace Nektar
{
    namespace MultiRegions
    {
        // Forward declarations
        class GlobalLinSys;
        class AssemblyMapDG;

        class AssemblyMapCG;
        class GlobalLinSysKey;
        class GlobalMatrix;

        enum Direction
        {
            eX,
            eY,
            eZ,
            eS,
            eN
        };

        enum ExpansionType
        {
            e0D,
            e1D,
            e2D,
            e3DH1D,
            e3DH2D,
            e3D,
            eNoType
        };       
        
        MultiRegions::Direction const DirCartesianMap[] =
        {
            eX,
            eY,
            eZ
        }; 
    
        /// A map between global matrix keys and their associated block
        /// matrices.
        typedef std::map<GlobalMatrixKey,DNekScalBlkMatSharedPtr> BlockMatrixMap;
        /// A shared pointer to a BlockMatrixMap.
        typedef boost::shared_ptr<BlockMatrixMap> BlockMatrixMapShPtr;
                   

        /// Base class for all multi-elemental spectral/hp expansions.
        class ExpList: public boost::enable_shared_from_this<ExpList>
        {
        public:
            /// The default constructor.
            MULTI_REGIONS_EXPORT ExpList();
            
            /// The default constructor.
            MULTI_REGIONS_EXPORT ExpList(
                    const LibUtilities::SessionReaderSharedPtr &pSession);

            /// The default constructor.
            MULTI_REGIONS_EXPORT ExpList(
                    const LibUtilities::SessionReaderSharedPtr &pSession,
                    const SpatialDomains::MeshGraphSharedPtr &pGraph);
            
            /// Constructor copying only elements defined in eIds.
            MULTI_REGIONS_EXPORT ExpList(
                const ExpList &in,
                const std::vector<unsigned int> &eIDs,
                const bool DeclareCoeffPhysArrays = true);

            /// The copy constructor.
            MULTI_REGIONS_EXPORT ExpList(
                const ExpList &in,
                const bool DeclareCoeffPhysArrays = true);

            /// The default destructor.
            MULTI_REGIONS_EXPORT virtual ~ExpList();
            
            /// Returns the total number of local degrees of freedom
            /// \f$N_{\mathrm{eof}}=\sum_{e=1}^{{N_{\mathrm{el}}}}N^{e}_m\f$.
            inline int GetNcoeffs(void) const;
            
            /// Returns the total number of local degrees of freedom
            /// for element eid
            MULTI_REGIONS_EXPORT int GetNcoeffs(const int eid) const;

            /// Returns the type of the expansion
            MULTI_REGIONS_EXPORT ExpansionType GetExpType(void);
            
            /// Returns the type of the expansion
            MULTI_REGIONS_EXPORT void SetExpType(ExpansionType Type);
    
            /// Evaulates the maximum number of modes in the elemental basis
            /// order over all elements
            inline int EvalBasisNumModesMax(void) const;

            /// Returns the vector of the number of modes in the elemental
            /// basis order over all elements.
            MULTI_REGIONS_EXPORT const Array<OneD,int>
                EvalBasisNumModesMaxPerExp(void) const;

            /// Returns the total number of quadrature points #m_npoints
            /// \f$=Q_{\mathrm{tot}}\f$.
            inline int GetTotPoints(void) const;

            /// Returns the total number of quadrature points for eid's element
            /// \f$=Q_{\mathrm{tot}}\f$.
            inline int GetTotPoints(const int eid) const;
            
            /// Returns the total number of quadrature points #m_npoints
            /// \f$=Q_{\mathrm{tot}}\f$.
            inline int GetNpoints(void) const;

            
            /// Returns the total number of qudature points scaled by
            /// the factor scale on each 1D direction
            inline int Get1DScaledTotPoints(const NekDouble scale) const;
            
            /// Sets the wave space to the one of the possible configuration
            /// true or false
            inline void SetWaveSpace(const bool wavespace);
            
            
            ///Set Modified Basis for the stability analysis
            inline void SetModifiedBasis(const bool modbasis);
            
            /// Set the \a i th value of \a m_phys to value \a val
            inline void SetPhys(int i, NekDouble val);
            
            /// This function returns the third direction expansion condition,
            /// which can be in wave space (coefficient) or not
            /// It is stored in the variable m_WaveSpace.
            inline bool GetWaveSpace(void) const;
            
            /// Fills the array #m_phys
            inline void SetPhys(const Array<OneD, const NekDouble> &inarray);

            /// Sets the array #m_phys
            inline void SetPhysArray(Array<OneD, NekDouble> &inarray);

            /// This function manually sets whether the array of physical
            /// values \f$\boldsymbol{u}_l\f$ (implemented as #m_phys) is
            /// filled or not.
            inline void SetPhysState(const bool physState);

            /// This function indicates whether the array of physical values
            /// \f$\boldsymbol{u}_l\f$ (implemented as #m_phys) is filled or
            /// not.
            inline bool GetPhysState(void) const;

            /// This function integrates a function \f$f(\boldsymbol{x})\f$
            /// over the domain consisting of all the elements of the expansion.
            MULTI_REGIONS_EXPORT NekDouble PhysIntegral (void);

            /// This function integrates a function \f$f(\boldsymbol{x})\f$
            /// over the domain consisting of all the elements of the expansion.
            MULTI_REGIONS_EXPORT NekDouble PhysIntegral(
                const Array<OneD,
                const NekDouble> &inarray);

            /// This function calculates the inner product of a function
            /// \f$f(\boldsymbol{x})\f$ with respect to all \emph{local}
            /// expansion modes \f$\phi_n^e(\boldsymbol{x})\f$.
            inline void   IProductWRTBase_IterPerExp(
                const Array<OneD, const NekDouble> &inarray,
                      Array<OneD,       NekDouble> &outarray);

            ///
            inline void IProductWRTBase(
                const Array<OneD, const NekDouble> &inarray,
                      Array<OneD,       NekDouble> &outarray,
                      CoeffState coeffstate = eLocal);

            /// This function calculates the inner product of a function
            /// \f$f(\boldsymbol{x})\f$ with respect to the derivative (in
            /// direction \param dir) of all \emph{local} expansion modes
            /// \f$\phi_n^e(\boldsymbol{x})\f$.
            MULTI_REGIONS_EXPORT void   IProductWRTDerivBase
                (const int dir,
                 const Array<OneD, const NekDouble> &inarray,
                       Array<OneD,       NekDouble> &outarray);

            /// This function calculates the inner product of a function
            /// \f$f(\boldsymbol{x})\f$ with respect to the derivative (in
            /// direction \param dir) of all \emph{local} expansion modes
            /// \f$\phi_n^e(\boldsymbol{x})\f$.
            MULTI_REGIONS_EXPORT void   IProductWRTDerivBase
                (const Array<OneD, const Array<OneD, NekDouble> > &inarray,
                 Array<OneD,       NekDouble> &outarray);

            /// This function elementally evaluates the forward transformation
            /// of a function \f$u(\boldsymbol{x})\f$ onto the global
            /// spectral/hp expansion.
            inline void  FwdTrans_IterPerExp (
                const Array<OneD,
                const NekDouble> &inarray,
                      Array<OneD,NekDouble> &outarray);

            ///
            inline void FwdTrans(
                const Array<OneD,
                const NekDouble> &inarray,
                      Array<OneD,       NekDouble> &outarray,
                      CoeffState coeffstate = eLocal);

            /// This function elementally mulplies the coefficient space of
            /// Sin my the elemental inverse of the mass matrix.
            MULTI_REGIONS_EXPORT void  MultiplyByElmtInvMass (
                 const Array<OneD,
                 const NekDouble> &inarray,
                 Array<OneD,       NekDouble> &outarray);

            ///
            inline void MultiplyByInvMassMatrix(
                const Array<OneD,const NekDouble> &inarray,
                      Array<OneD,      NekDouble> &outarray,
                      CoeffState coeffstate = eLocal);

            /// Smooth a field across elements
            inline void SmoothField(Array<OneD,NekDouble> &field);

            /// Solve helmholtz problem
            inline void HelmSolve(
                const Array<OneD, const NekDouble> &inarray,
                      Array<OneD,       NekDouble> &outarray,
                const FlagList &flags,
                const StdRegions::ConstFactorMap &factors,
                const StdRegions::VarCoeffMap &varcoeff =
                                StdRegions::NullVarCoeffMap,
                const Array<OneD, const NekDouble> &dirForcing =
                NullNekDouble1DArray,
                const bool PhysSpaceForcing = true);
            
            /// Solve Advection Diffusion Reaction
            inline void LinearAdvectionDiffusionReactionSolve(
                const Array<OneD, Array<OneD, NekDouble> > &velocity,
                       const Array<OneD, const NekDouble> &inarray,
                       Array<OneD, NekDouble> &outarray,
                       const NekDouble lambda,
                       CoeffState coeffstate = eLocal, 
                       const Array<OneD, const NekDouble>&
                       dirForcing = NullNekDouble1DArray);


            /// Solve Advection Diffusion Reaction
            inline void LinearAdvectionReactionSolve(
                const Array<OneD, Array<OneD, NekDouble> > &velocity,
                const Array<OneD, const NekDouble> &inarray,
                      Array<OneD, NekDouble> &outarray,
                const NekDouble lambda,
                      CoeffState coeffstate = eLocal,
                const Array<OneD, const NekDouble>&
                      dirForcing = NullNekDouble1DArray);

            ///
            MULTI_REGIONS_EXPORT void FwdTrans_BndConstrained(
                const Array<OneD, const NekDouble> &inarray,
                      Array<OneD,       NekDouble> &outarray);


            /// This function elementally evaluates the backward transformation
            /// of the global spectral/hp element expansion.
            inline void BwdTrans_IterPerExp (
                const Array<OneD, const NekDouble> &inarray,
                      Array<OneD,NekDouble> &outarray);

            ///
            inline void BwdTrans (
                const Array<OneD,
                const NekDouble> &inarray,
                      Array<OneD,NekDouble> &outarray,
                      CoeffState coeffstate = eLocal);

            /// This function calculates the coordinates of all the elemental
            /// quadrature points \f$\boldsymbol{x}_i\f$.
            inline void GetCoords(
                Array<OneD, NekDouble> &coord_0,
                Array<OneD, NekDouble> &coord_1 = NullNekDouble1DArray,
                Array<OneD, NekDouble> &coord_2 = NullNekDouble1DArray);
            
            // Homogeneous transforms
            inline void HomogeneousFwdTrans(
                const Array<OneD, const NekDouble> &inarray,
                      Array<OneD, NekDouble> &outarray,
                      CoeffState coeffstate = eLocal,
                bool Shuff = true,
                bool UnShuff = true);
            
            inline void HomogeneousBwdTrans(
                const Array<OneD, const NekDouble> &inarray,
                      Array<OneD, NekDouble> &outarray,
                      CoeffState coeffstate = eLocal,
                bool Shuff = true,
                bool UnShuff = true);
            
            inline void DealiasedProd(
                const Array<OneD, NekDouble> &inarray1,
                const Array<OneD, NekDouble> &inarray2,
                      Array<OneD, NekDouble> &outarray,
                      CoeffState coeffstate = eLocal);

            inline void DealiasedDotProd(
                const Array<OneD, Array<OneD, NekDouble> > &inarray1,
                const Array<OneD, Array<OneD, NekDouble> > &inarray2,
                      Array<OneD, Array<OneD, NekDouble> > &outarray,
                      CoeffState coeffstate = eLocal);

            inline void GetBCValues(
                      Array<OneD, NekDouble> &BndVals,
                const Array<OneD, NekDouble> &TotField,
                int BndID);
            
            inline void NormVectorIProductWRTBase(
                Array<OneD, const NekDouble> &V1,
                Array<OneD, const NekDouble> &V2,
                Array<OneD, NekDouble> &outarray,
                int BndID);
            
            inline void NormVectorIProductWRTBase(
                Array<OneD, Array<OneD, NekDouble> > &V,
                Array<OneD, NekDouble> &outarray);
            
            /// Apply geometry information to each expansion.
            MULTI_REGIONS_EXPORT void ApplyGeomInfo();

            /// Reset geometry information and reset matrices
            MULTI_REGIONS_EXPORT void Reset()
            {
                v_Reset();
            }

            void WriteTecplotHeader(std::ostream &outfile,
                                    std::string var = "")
            {
                v_WriteTecplotHeader(outfile, var);
            }

            void WriteTecplotZone(
                std::ostream &outfile,
                int expansion = -1)
            {
                v_WriteTecplotZone(outfile, expansion);
            }

            void WriteTecplotField(std::ostream &outfile,
                                   int expansion = -1)
            {
                v_WriteTecplotField(outfile, expansion);
            }

            void WriteTecplotConnectivity(std::ostream &outfile,
                                          int expansion = -1)
            {
                v_WriteTecplotConnectivity(outfile, expansion);
            }

            MULTI_REGIONS_EXPORT void WriteVtkHeader(std::ostream &outfile);
            MULTI_REGIONS_EXPORT void WriteVtkFooter(std::ostream &outfile);

            void WriteVtkPieceHeader(std::ostream &outfile, int expansion,
                                     int istrip = 0)
            {
                v_WriteVtkPieceHeader(outfile, expansion, istrip);
            }

            MULTI_REGIONS_EXPORT void WriteVtkPieceFooter(
                std::ostream &outfile,
                int expansion);

            void WriteVtkPieceData  (
                std::ostream &outfile,
                int expansion,
                std::string var = "v")
            {
                v_WriteVtkPieceData(outfile, expansion, var);
            }

            /// This function returns the dimension of the coordinates of the
            /// element \a eid.
            // inline
            MULTI_REGIONS_EXPORT int GetCoordim(int eid);

            /// Set the \a i th coefficiient in \a m_coeffs to value \a val
            inline void SetCoeff(int i, NekDouble val);
            
            /// Set the \a i th coefficiient in  #m_coeffs to value \a val
            inline void SetCoeffs(int i, NekDouble val);

            /// Set the  #m_coeffs array to inarray
            inline void SetCoeffsArray(Array<OneD, NekDouble> &inarray);

            /// This function returns (a reference to) the array
            /// \f$\boldsymbol{\hat{u}}_l\f$ (implemented as #m_coeffs)
            /// containing all local expansion coefficients.
            inline const Array<OneD, const NekDouble> &GetCoeffs() const;

            /// Impose Dirichlet Boundary Conditions onto Array
            inline void ImposeDirichletConditions(
                Array<OneD,NekDouble>& outarray);


            /// Fill Bnd Condition expansion from the values stored in expansion
            inline void FillBndCondFromField(void);

            /// Fill Bnd Condition expansion in nreg from the values stored in expansion
            inline void FillBndCondFromField(const int nreg);

            /// Gathers the global coefficients \f$\boldsymbol{\hat{u}}_g\f$
            /// from the local coefficients \f$\boldsymbol{\hat{u}}_l\f$.
            // inline
            MULTI_REGIONS_EXPORT inline void LocalToGlobal(bool useComm = true);

            MULTI_REGIONS_EXPORT inline void LocalToGlobal(
                const Array<OneD, const NekDouble> &inarray,
                Array<OneD,NekDouble> &outarray,
                bool useComm = true);

            /// Scatters from the global coefficients
            /// \f$\boldsymbol{\hat{u}}_g\f$ to the local coefficients
            /// \f$\boldsymbol{\hat{u}}_l\f$.
            // inline
            MULTI_REGIONS_EXPORT inline void GlobalToLocal(void);

            /**
             * This operation is evaluated as:
             * \f{tabbing}
             * \hspace{1cm}  \= Do \= $e=$  $1, N_{\mathrm{el}}$ \\
             * \> \> Do \= $i=$  $0,N_m^e-1$ \\
             * \> \> \> $\boldsymbol{\hat{u}}^{e}[i] = \mbox{sign}[e][i] \cdot
             * \boldsymbol{\hat{u}}_g[\mbox{map}[e][i]]$ \\
             * \> \> continue \\
             * \> continue
             * \f}
             * where \a map\f$[e][i]\f$ is the mapping array and \a
             * sign\f$[e][i]\f$ is an array of similar dimensions ensuring the
             * correct modal connectivity between the different elements (both
             * these arrays are contained in the data member #m_locToGloMap). This
             * operation is equivalent to the scatter operation
             * \f$\boldsymbol{\hat{u}}_l=\mathcal{A}\boldsymbol{\hat{u}}_g\f$,
             * where \f$\mathcal{A}\f$ is the
             * \f$N_{\mathrm{eof}}\times N_{\mathrm{dof}}\f$ permutation matrix.
             *
             * @param   inarray     An array of size \f$N_\mathrm{dof}\f$
             *                      containing the global degrees of freedom
             *                      \f$\boldsymbol{x}_g\f$.
             * @param   outarray    The resulting local degrees of freedom
             *                      \f$\boldsymbol{x}_l\f$ will be stored in this
             *                      array of size \f$N_\mathrm{eof}\f$.
             */
            MULTI_REGIONS_EXPORT inline void GlobalToLocal(
                const Array<OneD, const NekDouble> &inarray,
                Array<OneD,NekDouble> &outarray);

            /// Get the \a i th value  (coefficient) of #m_coeffs
            inline NekDouble GetCoeff(int i);

            /// Get the \a i th value  (coefficient) of #m_coeffs
            inline NekDouble GetCoeffs(int i);

            /// This function returns (a reference to) the array
            /// \f$\boldsymbol{u}_l\f$ (implemented as #m_phys) containing the
            /// function \f$u^{\delta}(\boldsymbol{x})\f$ evaluated at the
            /// quadrature points.
            // inline
            MULTI_REGIONS_EXPORT const Array<OneD, const NekDouble>
                &GetPhys()  const;

            /// This function calculates the \f$L_\infty\f$ error of the global
            /// spectral/hp element approximation.
            MULTI_REGIONS_EXPORT NekDouble Linf (
                const Array<OneD, const NekDouble> &inarray,
                const Array<OneD, const NekDouble> &soln = NullNekDouble1DArray);

            /// This function calculates the \f$L_2\f$ error with
            /// respect to soln of the global
            /// spectral/hp element approximation.
            NekDouble L2(
                const Array<OneD, const NekDouble> &inarray,
                const Array<OneD, const NekDouble> &soln = NullNekDouble1DArray)
            {
                return v_L2(inarray, soln);
            }

            /// Calculates the \f$H^1\f$ error of the global spectral/hp
            /// element approximation.
            MULTI_REGIONS_EXPORT NekDouble H1 (
                const Array<OneD, const NekDouble> &inarray,
                const Array<OneD, const NekDouble> &soln = NullNekDouble1DArray);
            
            NekDouble Integral (const Array<OneD, const NekDouble> &inarray)
            {
                return v_Integral(inarray);
            }

            /// This function calculates the energy associated with
            /// each one of the modesof a 3D homogeneous nD expansion
            Array<OneD, const NekDouble> HomogeneousEnergy (void)
            {
                return v_HomogeneousEnergy();
            }

            /// This function sets the Spectral Vanishing Viscosity 
            /// in homogeneous1D expansion. 
            void SetHomo1DSpecVanVisc(Array<OneD, NekDouble> visc)
            {
                v_SetHomo1DSpecVanVisc(visc);
            }

            /// This function returns a vector containing the wave
            /// numbers in z-direction associated
            /// with the 3D homogenous expansion. Required if a
            /// parellelisation is applied in the Fourier direction
            Array<OneD, const unsigned int> GetZIDs(void)
            {
                return v_GetZIDs();
            }

            /// This function returns the transposition class
            /// associaed with the homogeneous expansion.
            LibUtilities::TranspositionSharedPtr GetTransposition(void)
            {
                return v_GetTransposition();
            }
            
            /// This function returns the Width of homogeneous direction
            /// associaed with the homogeneous expansion.
            NekDouble GetHomoLen(void)
            {
                return v_GetHomoLen();
            }
            
            /// This function returns a vector containing the wave
            /// numbers in y-direction associated
            /// with the 3D homogenous expansion. Required if a
            /// parellelisation is applied in the Fourier direction
            Array<OneD, const unsigned int> GetYIDs(void)
            {
                return v_GetYIDs();
            }

            /// This function interpolates the physical space points in
            /// \a inarray to \a outarray using the same points defined in the
            /// expansion but where the number of points are rescaled
            /// by \a 1DScale
            void PhysInterp1DScaled(
                const NekDouble scale,
                const Array<OneD, NekDouble> &inarray,
                      Array<OneD, NekDouble>  &outarray)
            {
                v_PhysInterp1DScaled(scale, inarray,outarray);                
            }

            /// This function Galerkin projects the physical space points in
            /// \a inarray to \a outarray where inarray is assumed to
            /// be defined in the expansion but where the number of
            /// points are rescaled by \a 1DScale
            void PhysGalerkinProjection1DScaled(
                const NekDouble scale,
                const Array<OneD, NekDouble> &inarray,
                      Array<OneD, NekDouble> &outarray)
            {
                v_PhysGalerkinProjection1DScaled(scale, inarray, outarray);
            } 

            /// This function returns the number of elements in the expansion.
            inline int GetExpSize(void);


            /// This function returns the number of elements in the
            /// expansion which may be different for a homogeoenous extended
            /// expansionp.
            inline int GetNumElmts(void)
            {
                return v_GetNumElmts();
            }

            /// This function returns the vector of elements in the expansion.
            inline const boost::shared_ptr<LocalRegions::ExpansionVector>
                    GetExp() const;

            /// This function returns (a shared pointer to) the local elemental
            /// expansion of the \f$n^{\mathrm{th}}\f$ element.
            inline LocalRegions::ExpansionSharedPtr& GetExp(int n) const;

            /// This function returns (a shared pointer to) the local elemental
            /// expansion containing the arbitrary point given by \a gloCoord.
            MULTI_REGIONS_EXPORT LocalRegions::ExpansionSharedPtr& GetExp(
                const Array<OneD, const NekDouble> &gloCoord);

            /** This function returns the index of the local elemental
             * expansion containing the arbitrary point given by \a gloCoord.
             **/
            MULTI_REGIONS_EXPORT int GetExpIndex(
                const Array<OneD, const NekDouble> &gloCoord,
                NekDouble tol = 0.0,
                bool returnNearestElmt = false);

            /** This function returns the index and the Local
             * Cartesian Coordinates \a locCoords of the local
             * elemental expansion containing the arbitrary point
             * given by \a gloCoords.
             **/ 
            MULTI_REGIONS_EXPORT int GetExpIndex(
                const Array<OneD, const NekDouble> &gloCoords, 
                Array<OneD, NekDouble>       &locCoords,
                NekDouble tol = 0.0,
                bool returnNearestElmt = false);

            /// Get the start offset position for a global list of #m_coeffs
            /// correspoinding to element n.
            inline int GetCoeff_Offset(int n) const;

            /// Get the start offset position for a global list of m_phys
            /// correspoinding to element n.
            inline int GetPhys_Offset(int n) const;

            /// Get the element id associated with the n th
            /// consecutive block of data in  #m_phys and #m_coeffs
            inline int GetOffset_Elmt_Id(int n) const;

            /// This function returns (a reference to) the array
            /// \f$\boldsymbol{\hat{u}}_l\f$ (implemented as #m_coeffs)
            /// containing all local expansion coefficients.
            inline Array<OneD, NekDouble> &UpdateCoeffs();

            /// This function returns (a reference to) the array
            /// \f$\boldsymbol{u}_l\f$ (implemented as #m_phys) containing the
            /// function \f$u^{\delta}(\boldsymbol{x})\f$ evaluated at the
            /// quadrature points.
            inline Array<OneD, NekDouble> &UpdatePhys();

            inline void PhysDeriv(
                Direction edir,
                const Array<OneD, const NekDouble> &inarray,
                      Array<OneD, NekDouble> &out_d);    
        
            /// This function discretely evaluates the derivative of a function
            /// \f$f(\boldsymbol{x})\f$ on the domain consisting of all
            /// elements of the expansion.
            inline 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);
            
            inline void PhysDeriv(
                const int dir,
                const Array<OneD, const NekDouble> &inarray,
                      Array<OneD, NekDouble> &out_d);

            inline void CurlCurl(
                Array<OneD, Array<OneD, NekDouble> > &Vel,
                Array<OneD, Array<OneD, NekDouble> > &Q);

            // functions associated with DisContField
            inline const Array<OneD, const  boost::shared_ptr<ExpList> >
                &GetBndCondExpansions();
            
            inline boost::shared_ptr<ExpList> &UpdateBndCondExpansion(int i);
            
            inline void Upwind(
                const Array<OneD, const Array<OneD,       NekDouble> > &Vec,
                const Array<OneD,                   const NekDouble>   &Fwd,
                const Array<OneD,                   const NekDouble>   &Bwd,
                      Array<OneD,                         NekDouble>   &Upwind);

            inline void Upwind(
                const Array<OneD, const NekDouble> &Vn,
                const Array<OneD, const NekDouble> &Fwd,
                const Array<OneD, const NekDouble> &Bwd,
                      Array<OneD,       NekDouble> &Upwind);
            
            /**
             * Return a reference to the trace space associated with this
             * expansion list.
             */
            inline boost::shared_ptr<ExpList> &GetTrace();
            
            inline boost::shared_ptr<AssemblyMapDG> &GetTraceMap(void);
            
            inline const Array<OneD, const int> &GetTraceBndMap(void);

            inline void GetNormals(Array<OneD, Array<OneD, NekDouble> > &normals);

            inline void AddTraceIntegral(
                const Array<OneD, const NekDouble> &Fx,
                const Array<OneD, const NekDouble> &Fy,
                Array<OneD, NekDouble> &outarray);

            inline void AddTraceIntegral(
                const Array<OneD, const NekDouble> &Fn,
                      Array<OneD, NekDouble> &outarray);

            inline void AddFwdBwdTraceIntegral(
                const Array<OneD, const NekDouble> &Fwd,
                const Array<OneD, const NekDouble> &Bwd,
                      Array<OneD, NekDouble> &outarray);

            inline void GetFwdBwdTracePhys(
                Array<OneD,NekDouble> &Fwd,
                Array<OneD,NekDouble> &Bwd);

            inline void GetFwdBwdTracePhys(
                const Array<OneD,const NekDouble> &field,
                      Array<OneD,NekDouble> &Fwd,
                      Array<OneD,NekDouble> &Bwd);

            inline const std::vector<bool> &GetLeftAdjacentFaces(void) const;
            
            inline void ExtractTracePhys(Array<OneD,NekDouble> &outarray);

            inline void ExtractTracePhys(
                const Array<OneD, const NekDouble> &inarray,
                      Array<OneD,NekDouble> &outarray);

            inline const Array<OneD, const SpatialDomains::
                BoundaryConditionShPtr>& GetBndConditions();

            inline Array<OneD, SpatialDomains::
                BoundaryConditionShPtr>& UpdateBndConditions();

            inline void EvaluateBoundaryConditions(
                const NekDouble   time      = 0.0,
                const std::string varName   = "",
                const             NekDouble = NekConstants::kNekUnsetDouble,
                const             NekDouble = NekConstants::kNekUnsetDouble);

            // Routines for continous matrix solution
            /// This function calculates the result of the multiplication of a
            /// matrix of type specified by \a mkey with a vector given by \a
            /// inarray.
            inline void GeneralMatrixOp(
                const GlobalMatrixKey             &gkey,
                const Array<OneD,const NekDouble> &inarray,
                      Array<OneD,      NekDouble> &outarray,
                      CoeffState coeffstate = eLocal);

            MULTI_REGIONS_EXPORT void GeneralMatrixOp_IterPerExp(
                const GlobalMatrixKey      &gkey,
                const Array<OneD,const NekDouble> &inarray,
                      Array<OneD,      NekDouble> &outarray);

            inline void SetUpPhysNormals();

            inline void GetBoundaryToElmtMap(Array<OneD, int> &ElmtID,
                                             Array<OneD,int> &EdgeID);
            
            inline void GetBndElmtExpansion(int i,
                            boost::shared_ptr<ExpList> &result,
                            const bool DeclareCoeffPhysArrays = true);
            
            inline void ExtractElmtToBndPhys(int i,
                            Array<OneD, NekDouble> &elmt,
                            Array<OneD, NekDouble> &boundary);
            
            inline void ExtractPhysToBndElmt(int i,
                            const Array<OneD, const NekDouble> &phys,
                            Array<OneD, NekDouble> &bndElmt);

            inline void ExtractPhysToBnd(int i,
                            const Array<OneD, const NekDouble> &phys,
                            Array<OneD, NekDouble> &bnd);
            
            inline void GetBoundaryNormals(int i,
                            Array<OneD, Array<OneD, NekDouble> > &normals);

            MULTI_REGIONS_EXPORT void  GeneralGetFieldDefinitions(
                std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef,
                int NumHomoDir = 0,
                Array<OneD, LibUtilities::BasisSharedPtr> &HomoBasis =
                    LibUtilities::NullBasisSharedPtr1DArray,
                std::vector<NekDouble> &HomoLen =
                    LibUtilities::NullNekDoubleVector,
                bool  homoStrips = false,
                std::vector<unsigned int> &HomoSIDs =
                    LibUtilities::NullUnsignedIntVector,
                std::vector<unsigned int> &HomoZIDs =
                    LibUtilities::NullUnsignedIntVector,
                std::vector<unsigned int> &HomoYIDs =
                    LibUtilities::NullUnsignedIntVector);
            

            const NekOptimize::GlobalOptParamSharedPtr &GetGlobalOptParam(void)
            {
                return m_globalOptParam;
            }

            std::map<int, RobinBCInfoSharedPtr> GetRobinBCInfo()
            {
                return v_GetRobinBCInfo();
            }

            void GetPeriodicEntities(
                PeriodicMap &periodicVerts,
                PeriodicMap &periodicEdges,
                PeriodicMap &periodicFaces = NullPeriodicMap)
            {
                v_GetPeriodicEntities(periodicVerts, periodicEdges, periodicFaces);
            }

            std::vector<LibUtilities::FieldDefinitionsSharedPtr>
                GetFieldDefinitions()
            {
                return v_GetFieldDefinitions();
            }


            void GetFieldDefinitions(std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
            {
                v_GetFieldDefinitions(fielddef);
            }



            /// Append the element data listed in elements
            /// fielddef->m_ElementIDs onto fielddata
            void AppendFieldData(
                LibUtilities::FieldDefinitionsSharedPtr &fielddef,
                std::vector<NekDouble> &fielddata)
            {
                v_AppendFieldData(fielddef,fielddata);
            }

            
            /// Append the data in coeffs listed in elements
            /// fielddef->m_ElementIDs onto fielddata
            void AppendFieldData(
                LibUtilities::FieldDefinitionsSharedPtr &fielddef,
                std::vector<NekDouble> &fielddata,
                Array<OneD, NekDouble> &coeffs)
            {
                v_AppendFieldData(fielddef,fielddata,coeffs);
            }


            /** \brief Extract the data in fielddata into the coeffs
             * using the basic ExpList Elemental expansions rather
             * than planes in homogeneous case
             */ 
            MULTI_REGIONS_EXPORT void ExtractElmtDataToCoeffs(
                LibUtilities::FieldDefinitionsSharedPtr &fielddef,
                std::vector<NekDouble> &fielddata,
                std::string &field,
                Array<OneD, NekDouble> &coeffs);
            

            /** \brief Extract the data from fromField using
             * fromExpList the coeffs using the basic ExpList
             * Elemental expansions rather than planes in homogeneous
             * case
             */ 
            MULTI_REGIONS_EXPORT  void ExtractCoeffsToCoeffs(
                const boost::shared_ptr<ExpList> &fromExpList,
                const Array<OneD, const NekDouble> &fromCoeffs,
                      Array<OneD, NekDouble> &toCoeffs);
            
            
            //Extract data in fielddata into the m_coeffs_list for the 3D stability analysis (base flow is 2D)
            MULTI_REGIONS_EXPORT void ExtractDataToCoeffs(
                LibUtilities::FieldDefinitionsSharedPtr &fielddef,
                std::vector<NekDouble> &fielddata,
                std::string &field,
                Array<OneD, NekDouble> &coeffs);
            

            /// Returns a shared pointer to the current object.
            boost::shared_ptr<ExpList> GetSharedThisPtr()
            {
                return shared_from_this();
            }

            /// Returns the session object
            boost::shared_ptr<LibUtilities::SessionReader> GetSession() const
            {
                return m_session;
            }

            /// Returns the comm object
            boost::shared_ptr<LibUtilities::Comm> GetComm()
            {
                return m_comm;
            }

            SpatialDomains::MeshGraphSharedPtr GetGraph()
            {
                return m_graph;
            }

            // Wrapper functions for Homogeneous Expansions
            inline LibUtilities::BasisSharedPtr  GetHomogeneousBasis(void)
            {
                return v_GetHomogeneousBasis();
            }

            boost::shared_ptr<ExpList> &GetPlane(int n)
            {
                return v_GetPlane(n);
            }
           
            //expansion type
            ExpansionType m_expType;

            MULTI_REGIONS_EXPORT void CreateCollections(
                    Collections::ImplementationType ImpType
                                                    = Collections::eNoImpType);

            MULTI_REGIONS_EXPORT void ClearGlobalLinSysManager(void);

        protected:
            boost::shared_ptr<DNekMat> GenGlobalMatrixFull(
                const GlobalLinSysKey &mkey,
                const boost::shared_ptr<AssemblyMapCG> &locToGloMap);

            /// Communicator
            LibUtilities::CommSharedPtr m_comm;

            /// Session
            LibUtilities::SessionReaderSharedPtr m_session;

            /// Mesh associated with this expansion list.
            SpatialDomains::MeshGraphSharedPtr m_graph;

            /// The total number of local degrees of freedom. #m_ncoeffs
            /// \f$=N_{\mathrm{eof}}=\sum_{e=1}^{{N_{\mathrm{el}}}}N^{e}_l\f$
            int m_ncoeffs;

            /** The total number of quadrature points. #m_npoints
             *\f$=Q_{\mathrm{tot}}=\sum_{e=1}^{{N_{\mathrm{el}}}}N^{e}_Q\f$
             **/
            int m_npoints;
            
            /**
             * \brief Concatenation of all local expansion coefficients.
             *
             * The array of length #m_ncoeffs\f$=N_{\mathrm{eof}}\f$ which is
             * the concatenation of the local expansion coefficients
             * \f$\hat{u}_n^e\f$ over all \f$N_{\mathrm{el}}\f$ elements
             * \f[\mathrm{\texttt{m\_coeffs}}=\boldsymbol{\hat{u}}_{l} =
             * \underline{\boldsymbol{\hat{u}}}^e = \left [ \begin{array}{c}
             * \boldsymbol{\hat{u}}^{1} \       \
             * \boldsymbol{\hat{u}}^{2} \       \
             * \vdots \                                                 \
             * \boldsymbol{\hat{u}}^{{{N_{\mathrm{el}}}}} \end{array} \right ],
             * \quad
             * \mathrm{where}\quad \boldsymbol{\hat{u}}^{e}[n]=\hat{u}_n^{e}\f]
             */
            Array<OneD, NekDouble> m_coeffs;

            /**