EquationSystem.cpp

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////////
//
// File EquationSystem.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// 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: Main wrapper class for Advection Diffusion Reaction Solver
//
///////////////////////////////////////////////////////////////////////////////


#include <SolverUtils/EquationSystem.h>

#include <LocalRegions/MatrixKey.h>
#include <LibUtilities/BasicUtils/Equation.h>
#include <MultiRegions/ContField1D.h>
#include <MultiRegions/ContField2D.h>
#include <MultiRegions/ContField3D.h>
#include <MultiRegions/ContField3DHomogeneous1D.h>
#include <MultiRegions/ContField3DHomogeneous2D.h>

#include <MultiRegions/ExpList2D.h>     // for ExpList2D, etc
#include <MultiRegions/ExpList3D.h>     // for ExpList3D
#include <MultiRegions/ExpList3DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous2D.h>

#include <SolverUtils/AdvectionSystem.h>
#include <SolverUtils/Diffusion/Diffusion.h>

#include <GlobalMapping/Mapping.h>

#include <boost/format.hpp>
# include <boost/function.hpp>

#include <iostream>
#include <string>

using std::string;

namespace Nektar
{
    namespace SolverUtils
    {
        /**
         * @class EquationSystem
         *
         * This class is a base class for all solver implementations. It
         * provides the underlying generic functionality and interface for
         * solving equations.
         *
         * To solve a steady-state equation, create a derived class from this
         * class and reimplement the virtual functions to provide custom
         * implementation for the problem.
         *
         * To solve unsteady problems, derive from the UnsteadySystem class
         * instead which provides general time integration.
         */
        EquationSystemFactory& GetEquationSystemFactory()
        {
            typedef Loki::SingletonHolder<EquationSystemFactory,
                                          Loki::CreateUsingNew,
                                          Loki::NoDestroy,
                                          Loki::ClassLevelLockable> Type;
            return Type::Instance();
        }

        /**
         * This constructor is protected as the objects of this class are never
         * instantiated directly.
         * @param   pSession The session reader holding problem parameters.
         */
        EquationSystem::EquationSystem(
            const LibUtilities::SessionReaderSharedPtr& pSession)
            : m_comm (pSession->GetComm()),
              m_session (pSession),
              m_lambda (0),
              m_fieldMetaDataMap(LibUtilities::NullFieldMetaDataMap)
        {
            // set up session names in fieldMetaDataMap
            const vector<std::string> filenames = m_session->GetFilenames();
            
            for(int i = 0; i < filenames.size(); ++i)
            {
                string sessionname = "SessionName";
                sessionname += boost::lexical_cast<std::string>(i);
                m_fieldMetaDataMap[sessionname] = filenames[i];
            }
            
        }
        
        /**
         * @brief Initialisation object for EquationSystem.
         */
        void EquationSystem::v_InitObject()
        {
            // Save the basename of input file name for output details
            m_sessionName = m_session->GetSessionName();

            // Instantiate a field reader/writer
            m_fld = MemoryManager<LibUtilities::FieldIO>
                ::AllocateSharedPtr(
                    m_session->GetComm(),
                    m_session->DefinesCmdLineArgument("shared-filesystem"));

            // Read the geometry and the expansion information
            m_graph = SpatialDomains::MeshGraph::Read(m_session);

            // Also read and store the boundary conditions
            m_boundaryConditions =
                MemoryManager<SpatialDomains::BoundaryConditions>::
                    AllocateSharedPtr(m_session, m_graph);

            // Set space dimension for use in class
            m_spacedim = m_graph->GetSpaceDimension();
 
            // Setting parameteres for homogenous problems
            m_HomoDirec             = 0;
            m_useFFT                = false;
            m_homogen_dealiasing    = false;
            m_singleMode            = false;
            m_halfMode              = false;
            m_multipleModes         = false;
            m_HomogeneousType       = eNotHomogeneous;

            if (m_session->DefinesSolverInfo("HOMOGENEOUS"))
            {
                std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
                m_spacedim          = 3;

                if ((HomoStr == "HOMOGENEOUS1D") || (HomoStr == "Homogeneous1D")
                    || (HomoStr == "1D") || (HomoStr == "Homo1D"))
                {
                    m_HomogeneousType = eHomogeneous1D;
                    m_session->LoadParameter("LZ", m_LhomZ);
                    m_HomoDirec       = 1;
                
                    if(m_session->DefinesSolverInfo("ModeType"))
                    {
                        m_session->MatchSolverInfo("ModeType", "SingleMode", 
                                                   m_singleMode, false);
                        m_session->MatchSolverInfo("ModeType", "HalfMode", 
                                                   m_halfMode, false);
                        m_session->MatchSolverInfo("ModeType", "MultipleModes", 
                                                   m_multipleModes, false);
                    }

                    // Stability Analysis flags
                    if (m_session->DefinesSolverInfo("ModeType"))
                    {
                        if(m_singleMode)
                        {
                            m_npointsZ = 2;
                        }
                        else if(m_halfMode)
                        {
                            m_npointsZ = 1;
                        }
                        else if(m_multipleModes)
                        {
                            m_npointsZ = m_session->GetParameter("HomModesZ");
                        }
                        else
                        {
                            ASSERTL0(false, "SolverInfo ModeType not valid");
                        }
                    }
                    else 
                    {
                        m_npointsZ = m_session->GetParameter("HomModesZ");
                    }
                }

                if ((HomoStr == "HOMOGENEOUS2D") || (HomoStr == "Homogeneous2D")
                    || (HomoStr == "2D") || (HomoStr == "Homo2D"))
                {
                    m_HomogeneousType = eHomogeneous2D;
                    m_session->LoadParameter("HomModesY",   m_npointsY);
                    m_session->LoadParameter("LY",          m_LhomY);
                    m_session->LoadParameter("HomModesZ",   m_npointsZ);
                    m_session->LoadParameter("LZ",          m_LhomZ);
                    m_HomoDirec       = 2;
                }

                if ((HomoStr == "HOMOGENEOUS3D") || (HomoStr == "Homogeneous3D")
                    || (HomoStr == "3D") || (HomoStr == "Homo3D"))
                {
                    m_HomogeneousType = eHomogeneous3D;
                    m_session->LoadParameter("HomModesY",   m_npointsY);
                    m_session->LoadParameter("LY",          m_LhomY);
                    m_session->LoadParameter("HomModesZ",   m_npointsZ);
                    m_session->LoadParameter("LZ",          m_LhomZ);
                    m_HomoDirec       = 2;
                }

                m_session->MatchSolverInfo("USEFFT", "FFTW", m_useFFT, false);
            
                m_session->MatchSolverInfo("DEALIASING", "True", 
                                           m_homogen_dealiasing, false);
                if(m_homogen_dealiasing == false)
                {
                    m_session->MatchSolverInfo("DEALIASING", "On", 
                                               m_homogen_dealiasing, false);
                }
            }
            else
            {
                // set to default value so can use to identify 2d or 3D 
                // (homogeneous) expansions
                m_npointsZ = 1; 
            }
           
            m_session->MatchSolverInfo("SPECTRALHPDEALIASING", "True", 
                                       m_specHP_dealiasing, false);
            if (m_specHP_dealiasing == false)
            {
                m_session->MatchSolverInfo("SPECTRALHPDEALIASING", "On", 
                                           m_specHP_dealiasing, false);
            }
 
            // Options to determine type of projection from file or directly 
            // from constructor
            if (m_session->DefinesSolverInfo("PROJECTION"))
            {
                std::string ProjectStr = m_session->GetSolverInfo("PROJECTION");

                if ((ProjectStr == "Continuous") || (ProjectStr == "Galerkin") ||
                   (ProjectStr == "CONTINUOUS") || (ProjectStr == "GALERKIN"))
                {
                    m_projectionType = MultiRegions::eGalerkin;
                }
                else if ((ProjectStr == "MixedCGDG") ||
                        (ProjectStr == "Mixed_CG_Discontinuous"))
                {
                    m_projectionType = MultiRegions::eMixed_CG_Discontinuous;
                }                        
                else if(ProjectStr == "DisContinuous")
                {
                    m_projectionType = MultiRegions::eDiscontinuous;
                }
                else
                {
                    ASSERTL0(false,"PROJECTION value not recognised");
                }
            }
            else
            {
                cerr << "Projection type not specified in SOLVERINFO,"
                    "defaulting to continuous Galerkin" << endl;
                m_projectionType = MultiRegions::eGalerkin;
            }

            // Enforce singularity check for some problems
            m_checkIfSystemSingular = v_GetSystemSingularChecks();

            int i;
            int nvariables = m_session->GetVariables().size();
            bool DeclareCoeffPhysArrays = true;


            m_fields   = Array<OneD, MultiRegions::ExpListSharedPtr>(nvariables);
            m_spacedim = m_graph->GetSpaceDimension()+m_HomoDirec;
            m_expdim   = m_graph->GetMeshDimension();

            /// Continuous field
            if (m_projectionType == MultiRegions::eGalerkin ||
               m_projectionType == MultiRegions::eMixed_CG_Discontinuous)
            {
                switch(m_expdim)
                {
                case 1:
                    {
                        if (m_HomogeneousType == eHomogeneous2D
                            || m_HomogeneousType == eHomogeneous3D)
                        {
                            const LibUtilities::PointsKey PkeyY
                                (m_npointsY, LibUtilities::eFourierEvenlySpaced);
                            const LibUtilities::BasisKey  BkeyY
                                (LibUtilities::eFourier, m_npointsY, PkeyY);
                            const LibUtilities::PointsKey PkeyZ
                                (m_npointsZ, LibUtilities::eFourierEvenlySpaced);
                            const LibUtilities::BasisKey
                                BkeyZ(LibUtilities::eFourier, m_npointsZ, PkeyZ);

                            for (i = 0; i < m_fields.num_elements(); i++)
                            {
                                m_fields[i] = MemoryManager<MultiRegions
                                    ::ContField3DHomogeneous2D>
                                        ::AllocateSharedPtr(
                                            m_session, BkeyY, BkeyZ, m_LhomY, 
                                            m_LhomZ, m_useFFT, 
                                            m_homogen_dealiasing, m_graph, 
                                            m_session->GetVariable(i));
                            }
                        }
                        else
                        {
                            for (i = 0; i < m_fields.num_elements(); i++)
                            {
                                m_fields[i] = MemoryManager
                                <MultiRegions::ContField1D>::
                                    AllocateSharedPtr(
                                        m_session, m_graph,
                                        m_session->GetVariable(i));
                            }
                        }
                        break;
                    }
                case 2:
                    {
                        if (m_HomogeneousType == eHomogeneous1D)
                        {
                            // Fourier single mode stability analysis
                            if (m_singleMode)
                            {   
                                const LibUtilities::PointsKey PkeyZ(
                                    m_npointsZ,
                                    LibUtilities::eFourierSingleModeSpaced);
                                
                                const LibUtilities::BasisKey  BkeyZ(
                                    LibUtilities::eFourierSingleMode,
                                    m_npointsZ,
                                    PkeyZ);
                            
                                for(i = 0; i < m_fields.num_elements(); i++)
                                {
                                    m_fields[i] = MemoryManager<MultiRegions
                                        ::ContField3DHomogeneous1D>
                                            ::AllocateSharedPtr(
                                                m_session, BkeyZ, m_LhomZ, 
                                                m_useFFT, m_homogen_dealiasing, 
                                                m_graph, 
                                                m_session->GetVariable(i), 
                                                m_checkIfSystemSingular[i]);
                                }
                            }
                            // Half mode stability analysis
                            else if(m_halfMode)
                            {
                                const LibUtilities::PointsKey PkeyZ(
                                    m_npointsZ,
                                    LibUtilities::eFourierSingleModeSpaced);
                                    
                                const LibUtilities::BasisKey  BkeyZR(
                                    LibUtilities::eFourierHalfModeRe,
                                    m_npointsZ, PkeyZ);
                                
                                const LibUtilities::BasisKey  BkeyZI(
                                    LibUtilities::eFourierHalfModeIm,
                                    m_npointsZ, PkeyZ);
                                
                                    
                                for (i = 0; i < m_fields.num_elements(); i++)
                                {
                                    if(m_session->GetVariable(i).compare("w")
                                            == 0)
                                    {
                                        m_fields[i] = MemoryManager<MultiRegions
                                            ::ContField3DHomogeneous1D>
                                                ::AllocateSharedPtr(
                                                    m_session, BkeyZI, m_LhomZ, 
                                                    m_useFFT, 
                                                    m_homogen_dealiasing,
                                                    m_graph, 
                                                    m_session->GetVariable(i), 
                                                    m_checkIfSystemSingular[i]);
                                    }
                                    else
                                    {
                                        m_fields[i] = MemoryManager<MultiRegions
                                            ::ContField3DHomogeneous1D>
                                                ::AllocateSharedPtr(
                                                    m_session, BkeyZR, m_LhomZ, 
                                                    m_useFFT, m_homogen_dealiasing,
                                                    m_graph, 
                                                    m_session->GetVariable(i), 
                                                    m_checkIfSystemSingular[i]);
                                    }
                                    
                                    
                                }
                            }
                            // Normal homogeneous 1D
                            else
                            {   
                                const LibUtilities::PointsKey PkeyZ(
                                    m_npointsZ,
                                    LibUtilities::eFourierEvenlySpaced);
                                const LibUtilities::BasisKey  BkeyZ(
                                    LibUtilities::eFourier, m_npointsZ, PkeyZ);
                            
                                for (i = 0; i < m_fields.num_elements(); i++)
                                {
                                    m_fields[i] = MemoryManager<MultiRegions
                                        ::ContField3DHomogeneous1D>
                                            ::AllocateSharedPtr(
                                                m_session, BkeyZ, m_LhomZ, 
                                                m_useFFT, m_homogen_dealiasing,
                                                m_graph, 
                                                m_session->GetVariable(i), 
                                                m_checkIfSystemSingular[i]);
                                }
                            }
                        }
                        else
                        {
                            i = 0;
                            MultiRegions::ContField2DSharedPtr firstfield;
                            firstfield = MemoryManager<MultiRegions::
                                ContField2D>::AllocateSharedPtr(
                                    m_session, m_graph,
                                    m_session->GetVariable(i),
                                    DeclareCoeffPhysArrays,
                                    m_checkIfSystemSingular[0]);
                            m_fields[0] = firstfield;
                            for (i = 1; i < m_fields.num_elements(); i++)
                            {
                                if (m_graph->
                                      SameExpansions(m_session->GetVariable(0),
                                                     m_session->GetVariable(i)))
                                {
                                    m_fields[i] = MemoryManager<MultiRegions::
                                        ContField2D>::AllocateSharedPtr(
                                            *firstfield, m_graph,
                                            m_session->GetVariable(i),
                                            DeclareCoeffPhysArrays,
                                            m_checkIfSystemSingular[i]);
                                }
                                else
                                {
                                    m_fields[i] = MemoryManager<MultiRegions
                                        ::ContField2D>::AllocateSharedPtr(
                                            m_session, m_graph, 
                                            m_session->GetVariable(i),
                                            DeclareCoeffPhysArrays, 
                                            m_checkIfSystemSingular[i]);                                    
                                }
                            }

                            if (m_projectionType ==
                               MultiRegions::eMixed_CG_Discontinuous)
                            {
                                /// Setting up the normals
                                m_traceNormals =
                                    Array<OneD, Array<OneD, NekDouble> >
                                                                   (m_spacedim);
                                
                                for (i = 0; i < m_spacedim; ++i)
                                {
                                    m_traceNormals[i] = Array<OneD, NekDouble>
                                                            (GetTraceNpoints());
                                }
                                
                                m_fields[0]->GetTrace()->
                                    GetNormals(m_traceNormals);
                            }

                        }

                        break;
                    }
                case 3:
                    {
                        i = 0;
                        MultiRegions::ContField3DSharedPtr firstfield =
                            MemoryManager<MultiRegions::ContField3D>
                            ::AllocateSharedPtr(m_session, m_graph, 
                                                m_session->GetVariable(i),
                                                m_checkIfSystemSingular[i]);

                        m_fields[0] = firstfield;
                        for (i = 1; i < m_fields.num_elements(); i++)
                        {
                            if(m_graph->SameExpansions(
                                        m_session->GetVariable(0),
                                        m_session->GetVariable(i)))
                            {
                                m_fields[i] = MemoryManager<MultiRegions
                                    ::ContField3D>::AllocateSharedPtr(
                                        *firstfield, m_graph,
                                        m_session->GetVariable(i),
                                        m_checkIfSystemSingular[i]);
                            }
                            else
                            {
                                m_fields[i] = MemoryManager<MultiRegions
                                    ::ContField3D>::AllocateSharedPtr(
                                        m_session, m_graph, 
                                        m_session->GetVariable(i),
                                        m_checkIfSystemSingular[i]); 
                            }
                        }
                        
                        if (m_projectionType ==
                           MultiRegions::eMixed_CG_Discontinuous)
                        {
                            /// Setting up the normals
                            m_traceNormals =
                                Array<OneD, Array<OneD, NekDouble> >
                                                                (m_spacedim);
                            for(i = 0; i < m_spacedim; ++i)
                            {
                                m_traceNormals[i] =
                                    Array<OneD, NekDouble> (GetTraceNpoints());
                            }
                            
                            m_fields[0]->GetTrace()->GetNormals(m_traceNormals);
                            // Call the trace on all fields to ensure DG setup. 
                            for(i = 1; i < m_fields.num_elements(); ++i)
                            {
                                m_fields[i]->GetTrace();
                            }
                        }
                        break;
                    }
                default:
                    ASSERTL0(false,"Expansion dimension not recognised");
                    break;
                }
            }
            // Discontinuous field
            else
            {
                switch(m_expdim)
                {
                case 1:
                    {
                        if (m_HomogeneousType == eHomogeneous2D
                            || m_HomogeneousType == eHomogeneous3D)
                        {
                            const LibUtilities::PointsKey PkeyY(
                                m_npointsY, LibUtilities::eFourierEvenlySpaced);
                            const LibUtilities::BasisKey  BkeyY(
                                LibUtilities::eFourier, m_npointsY, PkeyY);
                            const LibUtilities::PointsKey PkeyZ(
                                m_npointsZ, LibUtilities::eFourierEvenlySpaced);
                            const LibUtilities::BasisKey  BkeyZ(
                                LibUtilities::eFourier, m_npointsZ, PkeyZ);

                            for (i = 0; i < m_fields.num_elements(); i++)
                            {
                                m_fields[i] = MemoryManager<MultiRegions
                                    ::DisContField3DHomogeneous2D>
                                        ::AllocateSharedPtr(
                                            m_session, BkeyY, BkeyZ, m_LhomY, 
                                            m_LhomZ, m_useFFT, 
                                            m_homogen_dealiasing, m_graph, 
                                            m_session->GetVariable(i));
                            }
                        }
                        else
                        {
                            for (i = 0; i < m_fields.num_elements(); i++)
                            {
                                m_fields[i] = MemoryManager<MultiRegions::
                                    DisContField1D>::AllocateSharedPtr(
                                        m_session, m_graph,
                                        m_session->GetVariable(i));
                            }
                        }

                        break;
                    }
                case 2:
                    {
                        if(m_HomogeneousType == eHomogeneous1D)
                        {
                            const LibUtilities::PointsKey PkeyZ(
                                m_npointsZ,LibUtilities::eFourierEvenlySpaced);
                            const LibUtilities::BasisKey BkeyZ(
                                LibUtilities::eFourier, m_npointsZ,PkeyZ);

                            for (i = 0; i < m_fields.num_elements(); i++)
                            {
                                m_fields[i] = MemoryManager<MultiRegions
                                    ::DisContField3DHomogeneous1D>
                                        ::AllocateSharedPtr(
                                            m_session, BkeyZ, m_LhomZ, m_useFFT, 
                                            m_homogen_dealiasing, m_graph, 
                                            m_session->GetVariable(i));
                            }
                        }
                        else
                        {
                            for (i = 0; i < m_fields.num_elements(); i++)
                            {
                                m_fields[i] = MemoryManager<MultiRegions::
                                    DisContField2D>::AllocateSharedPtr(
                                        m_session, m_graph,
                                        m_session->GetVariable(i));
                            }
                        }

                        break;
                    }
                case 3:
                    {
                        if (m_HomogeneousType == eHomogeneous3D)
                        {
                            ASSERTL0(false,
                              "3D fully periodic problems not implemented yet");
                        }
                        else
                        {
                            for (i = 0; i < m_fields.num_elements(); i++)
                            {
                                m_fields[i] = MemoryManager<MultiRegions::
                                    DisContField3D>::AllocateSharedPtr(
                                        m_session, m_graph,
                                        m_session->GetVariable(i));
                            }
                        }
                        break;
                    }
                    default:
                        ASSERTL0(false, "Expansion dimension not recognised");
                        break;
                }

                // Setting up the normals
                m_traceNormals =
                    Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
                
                for (i = 0; i < m_spacedim; ++i)
                {
                    m_traceNormals[i] =
                        Array<OneD, NekDouble> (GetTraceNpoints(), 0.0);
                }

                m_fields[0]->GetTrace()->GetNormals(m_traceNormals);
            }

            // Set Default Parameter
            m_session->LoadParameter("Time",          m_time,       0.0);
            m_session->LoadParameter("TimeStep",      m_timestep,   0.01);
            m_session->LoadParameter("NumSteps",      m_steps,      0);
            m_session->LoadParameter("IO_CheckSteps", m_checksteps, 0);
            m_session->LoadParameter("IO_CheckTime",  m_checktime,  0.0);
            m_session->LoadParameter("FinTime",       m_fintime,    0);
            m_session->LoadParameter("NumQuadPointsError",
                                     m_NumQuadPointsError, 0);

            m_nchk = 1;

            // Zero all physical fields initially
            ZeroPhysFields();
        }

        /**
         * @brief Destructor for class EquationSystem.
         */
        EquationSystem::~EquationSystem()
        {
            LibUtilities::NekManager<LocalRegions::MatrixKey,
                DNekScalMat, LocalRegions::MatrixKey::opLess>::ClearManager();
            LibUtilities::NekManager<LocalRegions::MatrixKey,
                DNekScalBlkMat, LocalRegions::MatrixKey::opLess>::ClearManager();
        }

        /**
         * Evaluates a physical function at each quadrature point in the domain.
         *
         * @param   pArray          The array into which to write the values.
         * @param   pEqn            The equation to evaluate.
         */
        void EquationSystem::EvaluateFunction(
            Array<OneD, Array<OneD, NekDouble> >& pArray,
            std::string pFunctionName,
            const NekDouble pTime,
            const int domain)
        {
            ASSERTL0(m_session->DefinesFunction(pFunctionName),
                     "Function '" + pFunctionName + "' does not exist.");

            std::vector<std::string> vFieldNames = m_session->GetVariables();

            for(int i = 0 ; i < vFieldNames.size(); i++)
            {
                EvaluateFunction(vFieldNames[i], pArray[i], pFunctionName,
                                 pTime, domain);
            }
        }

        /**
         * Populates a forcing function for each of the dependent variables 
         * using the expression provided by the BoundaryConditions object.
         * @param   force           Array of fields to assign forcing.
         */
        void EquationSystem::EvaluateFunction(
            std::vector<std::string> pFieldNames,
            Array<OneD, Array<OneD, NekDouble> > &pFields,
            const std::string& pFunctionName,
            const NekDouble& pTime,
            const int domain)
        {
            ASSERTL1(pFieldNames.size() == pFields.num_elements(),
                     "Function '" + pFunctionName
                     + "' variable list size mismatch with array storage.");
            ASSERTL0(m_session->DefinesFunction(pFunctionName),
                     "Function '" + pFunctionName + "' does not exist.");

            for(int i = 0; i < pFieldNames.size(); i++)
            {
                EvaluateFunction(pFieldNames[i], pFields[i], pFunctionName, pTime, domain);
            }
        }

        /**
         * Populates a function for each of the dependent variables using
         * the expression or filenames provided by the SessionReader object.
         * @param   force           Array of fields to assign forcing.
         */
        void EquationSystem::EvaluateFunction(
            std::vector<std::string> pFieldNames,
            Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
            const std::string& pFunctionName,
            const NekDouble& pTime,
            const int domain)
        {
            ASSERTL0(m_session->DefinesFunction(pFunctionName),
                     "Function '" + pFunctionName + "' does not exist.");
            ASSERTL0(pFieldNames.size() == pFields.num_elements(),
                     "Field list / name list size mismatch.");

            for(int i = 0; i < pFieldNames.size(); i++)
            {
                EvaluateFunction(pFieldNames[i], pFields[i]->UpdatePhys(),
                                 pFunctionName, pTime, domain);
                pFields[i]->FwdTrans_IterPerExp(pFields[i]->GetPhys(),
                                                pFields[i]->UpdateCoeffs());
            }

        }


        void EquationSystem::EvaluateFunction(
            std::string pFieldName,
            Array<OneD, NekDouble>& pArray,
            const std::string& pFunctionName,
            const NekDouble& pTime,
            const int domain)
        {
            ASSERTL0(m_session->DefinesFunction(pFunctionName),
                     "Function '" + pFunctionName + "' does not exist.");

            unsigned int nq = m_fields[0]->GetNpoints();
            if (pArray.num_elements() < nq)
            {
                pArray = Array<OneD, NekDouble>(nq);
            }

            LibUtilities::FunctionType vType;
            vType = m_session->GetFunctionType(pFunctionName, pFieldName,domain);
            if (vType == LibUtilities::eFunctionTypeExpression)
            {
                Array<OneD,NekDouble> x0(nq);
                Array<OneD,NekDouble> x1(nq);
                Array<OneD,NekDouble> x2(nq);
                
                // Get the coordinates (assuming all fields have the same
                // discretisation)
                m_fields[0]->GetCoords(x0,x1,x2);
                LibUtilities::EquationSharedPtr ffunc
                    = m_session->GetFunction(pFunctionName, pFieldName,domain);

                ffunc->Evaluate(x0,x1,x2,pTime,pArray);
            }
            else if (vType == LibUtilities::eFunctionTypeFile ||
                     vType == LibUtilities::eFunctionTypeTransientFile)
            {
                // check if we already read this pFunctionName + pFieldName
                // combination and stop processing if we are dealing with
                // a non-timedependent file
                std::string loadedKey = pFunctionName + pFieldName;
                if (m_loadedFields.count(loadedKey) != 0 && vType == LibUtilities::eFunctionTypeFile)
                {
                    return;
                }
                m_loadedFields.insert(loadedKey);

                std::string filename = m_session->GetFunctionFilename(
                    pFunctionName, pFieldName, domain);
                std::string fileVar = m_session->GetFunctionFilenameVariable(
                    pFunctionName, pFieldName, domain);

                if (fileVar.length() == 0)
                {
                    fileVar = pFieldName;
                }

                std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
                std::vector<std::vector<NekDouble> > FieldData;
                Array<OneD, NekDouble> vCoeffs(m_fields[0]->GetNcoeffs());
                Vmath::Zero(vCoeffs.num_elements(),vCoeffs,1);

                int numexp = m_fields[0]->GetExpSize();
                Array<OneD,int> ElementGIDs(numexp);

                // Define list of global element ids
                for(int i = 0; i < numexp; ++i)
                {
                    ElementGIDs[i] = m_fields[0]->GetExp(i)->GetGeom()->GetGlobalID();
                }

                //  In case of eFunctionTypeTransientFile, generate filename from
                //  format string
                if (vType == LibUtilities::eFunctionTypeTransientFile)
                {
                    try
                    {
#ifdef _WIN32
                        // We need this to make sure boost::format has always
                        // two digits in the exponents of Scientific notation.
                        unsigned int old_exponent_format;
                        old_exponent_format = _set_output_format(_TWO_DIGIT_EXPONENT);
                        filename = boost::str(boost::format(filename) % m_time);
                        _set_output_format(old_exponent_format);
#else
                        filename = boost::str(boost::format(filename) % m_time);
#endif
                    }
                    catch (...)
                    {
                        ASSERTL0(false, "Invalid Filename in function \""
                                + pFunctionName + "\", variable \"" + fileVar + "\"")
                    }
                }

                if (boost::filesystem::path(filename).extension() !=  ".pts")
                {
                    m_fld->Import(filename, FieldDef, FieldData,
                                LibUtilities::NullFieldMetaDataMap,
                                ElementGIDs);

                    int idx = -1;

                    // Loop over all the expansions
                    for (int i = 0; i < FieldDef.size(); ++i)
                    {
                        // Find the index of the required field in the
                        // expansion segment
                        for (int j = 0; j < FieldDef[i]->m_fields.size(); ++j)
                        {
                            if (FieldDef[i]->m_fields[j] == fileVar)
                            {
                                idx = j;
                            }
                        }

                        if (idx >= 0)
                        {
                            m_fields[0]->ExtractDataToCoeffs(
                                FieldDef[i], FieldData[i],
                                FieldDef[i]->m_fields[idx], vCoeffs);
                        }
                        else
                        {
                            cout << "Field " + fileVar + " not found." << endl;
                        }
                    }

                    m_fields[0]->BwdTrans_IterPerExp(vCoeffs, pArray);
                }
                else
                {
                    LibUtilities::PtsFieldSharedPtr ptsField;
                    LibUtilities::PtsIO ptsIO(m_session->GetComm());
                    ptsIO.Import(filename, ptsField);

                    Array <OneD,  Array<OneD,  NekDouble> > coords(3);
                    coords[0] = Array<OneD, NekDouble>(nq);
                    coords[1] = Array<OneD, NekDouble>(nq);
                    coords[2] = Array<OneD, NekDouble>(nq);
                    m_fields[0]->GetCoords(coords[0], coords[1], coords[2]);

                    //  check if we already computed this funcKey combination
                    std::string weightsKey = m_session->GetFunctionFilename(pFunctionName, pFieldName, domain);
                    if (m_interpWeights.count(weightsKey) != 0)
                    {
                        //  found, re-use
                        ptsField->SetWeights(m_interpWeights[weightsKey], m_interpInds[weightsKey]);
                    }
                    else
                    {
                        if (m_session->GetComm()->GetRank() == 0)
                        {
                            ptsField->setProgressCallback(&EquationSystem::PrintProgressbar, this);
                            cout << "Interpolating:       ";
                        }
                        ptsField->CalcWeights(coords);
                        if (m_session->GetComm()->GetRank() == 0)
                        {
                            cout << endl;
                        }
                        ptsField->GetWeights(m_interpWeights[weightsKey], m_interpInds[weightsKey]);
                    }

                    Array<OneD,  Array<OneD,  NekDouble> > intFields;
                    ptsField->Interpolate(intFields);

                    int fieldInd;
                    vector<string> fieldNames = ptsField->GetFieldNames();
                    for (fieldInd = 0; fieldInd < fieldNames.size(); ++fieldInd)
                    {
                        if (ptsField->GetFieldName(fieldInd) ==  pFieldName)
                        {
                            break;
                        }
                    }
                    ASSERTL0(fieldInd != fieldNames.size(),  "field not found");

                    pArray = intFields[fieldInd];
                }
            }
        }


        /**
         * @brief Provide a description of a function for a given field name.
         *
         * @param pFieldName     Field name.
         * @param pFunctionName  Function name.
         */
        std::string EquationSystem::DescribeFunction(
            std::string pFieldName,
            const std::string &pFunctionName,
            const int domain)
        {
            ASSERTL0(m_session->DefinesFunction(pFunctionName),
                     "Function '" + pFunctionName + "' does not exist.");
            
            std::string retVal;
            LibUtilities::FunctionType vType;
            
            vType = m_session->GetFunctionType(pFunctionName, pFieldName);
            if (vType == LibUtilities::eFunctionTypeExpression)
            {
                LibUtilities::EquationSharedPtr ffunc
                    = m_session->GetFunction(pFunctionName, pFieldName,domain);
                retVal = ffunc->GetExpression();
            }
            else if (vType == LibUtilities::eFunctionTypeFile)
            {
                std::string filename
                    = m_session->GetFunctionFilename(pFunctionName, pFieldName,domain);
                retVal = "from file " + filename;
            }
            
            return retVal;
        }
        
        /**
         * If boundary conditions are time-dependent, they will be evaluated at 
         * the time specified.
         * @param   time            The time at which to evaluate the BCs
         */
        void EquationSystem::SetBoundaryConditions(NekDouble time)
        {
            std::string varName;
            int nvariables = m_fields.num_elements();
            for (int i = 0; i < nvariables; ++i)
            {
                varName = m_session->GetVariable(i); 
                m_fields[i]->EvaluateBoundaryConditions(time, varName);
            }
        }

        /**
         * Compute the error in the L2-norm.
         * @param   field           The field to compare.
         * @param   exactsoln       The exact solution to compare with.
         * @param   Normalised      Normalise L2-error.
         * @returns                 Error in the L2-norm.
         */
        NekDouble EquationSystem::v_L2Error(
            unsigned int field,
            const Array<OneD, NekDouble> &exactsoln,
            bool Normalised)
        {       
            NekDouble L2error = -1.0;

            if (m_NumQuadPointsError == 0)
            {
                if (m_fields[field]->GetPhysState() == false)
                {
                    m_fields[field]->BwdTrans(m_fields[field]->GetCoeffs(),
                                              m_fields[field]->UpdatePhys());
                }

                if (exactsoln.num_elements())
                {
                    L2error = m_fields[field]->L2(m_fields[field]->GetPhys(), exactsoln);
                }
                else if (m_session->DefinesFunction("ExactSolution"))
                {
                    Array<OneD, NekDouble>
                        exactsoln(m_fields[field]->GetNpoints());

                    EvaluateFunction(m_session->GetVariable(field), exactsoln, 
                                     "ExactSolution", m_time);

                    L2error = m_fields[field]->L2(m_fields[field]->GetPhys(), exactsoln);
                }
                else
                {
                    L2error = m_fields[field]->L2(m_fields[field]->GetPhys());
                }

                if (Normalised == true)
                {
                    Array<OneD, NekDouble> one(m_fields[field]->GetNpoints(),
                                               1.0);

                    NekDouble Vol = m_fields[field]->PhysIntegral(one);
                    m_comm->AllReduce(Vol, LibUtilities::ReduceSum);

                    L2error = sqrt(L2error*L2error/Vol);
                }
            }
            else
            {
                Array<OneD,NekDouble> L2INF(2);
                L2INF = ErrorExtraPoints(field);
                L2error = L2INF[0];
            }
            return L2error;
        }

        /**
         * Compute the error in the L_inf-norm
         * @param   field           The field to compare.
         * @param   exactsoln       The exact solution to compare with.
         * @returns                 Error in the L_inft-norm.
         */
        NekDouble EquationSystem::v_LinfError (
                unsigned int field,
                const Array<OneD, NekDouble> &exactsoln)
        {
            NekDouble Linferror = -1.0;

            if (m_NumQuadPointsError == 0)
            {
                if (m_fields[field]->GetPhysState() == false)
                {
                    m_fields[field]->BwdTrans(m_fields[field]->GetCoeffs(),
                                              m_fields[field]->UpdatePhys());
                }

                if (exactsoln.num_elements())
                {
                    Linferror = m_fields[field]->Linf(m_fields[field]->GetPhys(), exactsoln);
                }
                else if (m_session->DefinesFunction("ExactSolution"))
                {
                    Array<OneD, NekDouble>
                        exactsoln(m_fields[field]->GetNpoints());

                    EvaluateFunction(m_session->GetVariable(field), exactsoln,
                                     "ExactSolution", m_time);

                    Linferror = m_fields[field]->Linf(m_fields[field]->GetPhys(), exactsoln);
                }
                else
                {
                    Linferror = m_fields[field]->Linf(m_fields[field]->GetPhys());
                }
            }
            else
            {
                Array<OneD,NekDouble> L2INF(2);
                L2INF = ErrorExtraPoints(field);
                Linferror = L2INF[1];
            }

            return Linferror;
        }

        /**
         * Compute the error in the L2-norm, L-inf for a larger number of 
         * quadrature points.
         * @param   field              The field to compare.
         * @returns                    Error in the L2-norm and L-inf norm.
         */
        Array<OneD,NekDouble> EquationSystem::ErrorExtraPoints(
            unsigned int field)
        {
            int NumModes = GetNumExpModes();
            Array<OneD,NekDouble> L2INF(2);

            const LibUtilities::PointsKey PkeyT1(
                m_NumQuadPointsError,LibUtilities::eGaussLobattoLegendre);
            const LibUtilities::PointsKey PkeyT2(
                m_NumQuadPointsError, LibUtilities::eGaussRadauMAlpha1Beta0);
            const LibUtilities::PointsKey PkeyQ1(
                m_NumQuadPointsError, LibUtilities::eGaussLobattoLegendre);
            const LibUtilities::PointsKey PkeyQ2(
                m_NumQuadPointsError, LibUtilities::eGaussLobattoLegendre);
            const LibUtilities::BasisKey  BkeyT1(
                LibUtilities::eModified_A,NumModes, PkeyT1);
            const LibUtilities::BasisKey  BkeyT2(
                LibUtilities::eModified_B, NumModes, PkeyT2);
            const LibUtilities::BasisKey  BkeyQ1(
                LibUtilities::eModified_A, NumModes, PkeyQ1);
            const LibUtilities::BasisKey  BkeyQ2(
                LibUtilities::eModified_A, NumModes, PkeyQ2);

            MultiRegions::ExpList2DSharedPtr ErrorExp =
                MemoryManager<MultiRegions::ExpList2D>::AllocateSharedPtr(
                        m_session, BkeyT1, BkeyT2, BkeyQ1, BkeyQ2, m_graph);

            int ErrorCoordim = ErrorExp->GetCoordim(0);
            int ErrorNq      = ErrorExp->GetTotPoints();

            Array<OneD,NekDouble> ErrorXc0(ErrorNq, 0.0);
            Array<OneD,NekDouble> ErrorXc1(ErrorNq, 0.0);
            Array<OneD,NekDouble> ErrorXc2(ErrorNq, 0.0);

            switch(ErrorCoordim)
            {
                case 1:
                    ErrorExp->GetCoords(ErrorXc0);
                    break;
                case 2:
                    ErrorExp->GetCoords(ErrorXc0, ErrorXc1);
                    break;
                case 3:
                    ErrorExp->GetCoords(ErrorXc0, ErrorXc1, ErrorXc2);
                    break;
            }
            LibUtilities::EquationSharedPtr exSol = 
                m_session->GetFunction("ExactSolution", field);
            
            // Evaluate the exact solution
            Array<OneD,NekDouble> ErrorSol(ErrorNq);

            exSol->Evaluate(ErrorXc0,ErrorXc1,ErrorXc2,m_time,ErrorSol);

            // Calcualte spectral/hp approximation on the quadrature points  
            // of this new expansion basis
            ErrorExp->BwdTrans_IterPerExp(m_fields[field]->GetCoeffs(), 
                                          ErrorExp->UpdatePhys());

            L2INF[0] = ErrorExp->L2  (ErrorExp->GetPhys(), ErrorSol);
            L2INF[1] = ErrorExp->Linf(ErrorExp->GetPhys(), ErrorSol);

            return L2INF;
        }


        /**
         * Set the physical fields based on a restart file, or a function
         * describing the initial condition given in the session.
         * @param  initialtime           Time at which to evaluate the function.
         * @param  dumpInitialConditions Write the initial condition to file?
         */
        void EquationSystem::v_SetInitialConditions(NekDouble initialtime,
                                                    bool dumpInitialConditions,
                                                    const int domain)
        {
            if (m_session->GetComm()->GetRank() == 0)
            {
                cout << "Initial Conditions:" << endl;
            }
        
            if (m_session->DefinesFunction("InitialConditions"))
            {
                EvaluateFunction(m_session->GetVariables(), m_fields, 
                                 "InitialConditions", m_time, domain);
                
                if (m_session->GetComm()->GetRank() == 0)
                {
                    
                    for (int i = 0; i < m_fields.num_elements(); ++i)
                    {
                        std::string varName = m_session->GetVariable(i);
                        cout << "  - Field " << varName << ": "
                             << DescribeFunction(varName, "InitialConditions",domain)
                             << endl;
                    }
                }
            }
            else
            {
                int nq = m_fields[0]->GetNpoints();
                for (int i = 0; i < m_fields.num_elements(); i++)
                {
                    Vmath::Zero(nq, m_fields[i]->UpdatePhys(), 1);
                    m_fields[i]->SetPhysState(true);
                    Vmath::Zero(m_fields[i]->GetNcoeffs(), 
                                m_fields[i]->UpdateCoeffs(), 1);
                    if (m_session->GetComm()->GetRank() == 0)
                    {
                        cout << "  - Field "    << m_session->GetVariable(i)
                             << ": 0 (default)" << endl;
                    }
                }

            }

            if (dumpInitialConditions && m_checksteps)
            {
                Checkpoint_Output(0);
            }
        }
    
    
        void EquationSystem::v_EvaluateExactSolution(
            unsigned int field,
            Array<OneD, NekDouble> &outfield,
            const NekDouble time)
        {
            ASSERTL0 (outfield.num_elements() == m_fields[field]->GetNpoints(),
                      "ExactSolution array size mismatch.");
            Vmath::Zero(outfield.num_elements(), outfield, 1);
            if (m_session->DefinesFunction("ExactSolution"))
            {
                EvaluateFunction(m_session->GetVariable(field), outfield, 
                                 "ExactSolution", time);
            }
        }


        /**
         * By default, nothing needs initialising at the EquationSystem level.
         */
        void EquationSystem::v_DoInitialise()
        {

        }
    
    
        void EquationSystem::InitialiseBaseFlow(
            Array<OneD, Array<OneD, NekDouble> > &base)
        {
            base = Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
            std::vector<std::string> vel;
            vel.push_back("Vx");
            vel.push_back("Vy");
            vel.push_back("Vz");
            vel.resize(m_spacedim);
            SetUpBaseFields(m_graph);
            EvaluateFunction(vel, base, "BaseFlow");
        }           
    
        void EquationSystem::SetUpBaseFields(
            SpatialDomains::MeshGraphSharedPtr &mesh)
        {
            int i;
            
            // The number of variables can be different from the dimension 
            // of the base flow
            int nvariables = m_session->GetVariables().size();
            m_base = Array<OneD, MultiRegions::ExpListSharedPtr>(nvariables);
            if (m_projectionType == MultiRegions::eGalerkin ||
                m_projectionType == MultiRegions::eMixed_CG_Discontinuous)
            {
                switch (m_expdim)
                {
                    case 1:
                    {
                        for(i = 0; i < m_base.num_elements(); i++)
                        {
                            m_base[i] = MemoryManager<MultiRegions::ContField1D>
                                ::AllocateSharedPtr(m_session, mesh, 
                                                    m_session->GetVariable(0));
                        }
                    }
                        break;
                    case 2:
                    {
                        if (m_HomogeneousType == eHomogeneous1D)
                        {
                            if (m_singleMode)
                            {
                                const LibUtilities::PointsKey PkeyZ(m_npointsZ,
                                        LibUtilities::eFourierSingleModeSpaced);
                                const LibUtilities::BasisKey  BkeyZ(
                                        LibUtilities::eFourier,
                                        m_npointsZ,PkeyZ);

                                for (i = 0 ; i < m_base.num_elements(); i++)
                                {
                                    m_base[i] = MemoryManager<MultiRegions
                                        ::ContField3DHomogeneous1D>
                                            ::AllocateSharedPtr(
                                                m_session, BkeyZ, m_LhomZ, 
                                                m_useFFT, m_homogen_dealiasing, 
                                                m_graph, 
                                                m_session->GetVariable(i));
                                    m_base[i]->SetWaveSpace(true);
                                }
                            }
                            else if (m_halfMode)
                            {
                                //1 plane field (half mode expansion)
                                const LibUtilities::PointsKey PkeyZ(m_npointsZ,
                                        LibUtilities::eFourierSingleModeSpaced);
                                const LibUtilities::BasisKey  BkeyZ(
                                        LibUtilities::eFourierHalfModeRe,
                                        m_npointsZ,PkeyZ);

                                for (i = 0 ; i < m_base.num_elements(); i++)
                                {
                                    m_base[i] = MemoryManager<MultiRegions
                                        ::ContField3DHomogeneous1D>
                                            ::AllocateSharedPtr(
                                                m_session, BkeyZ, m_LhomZ, 
                                                m_useFFT, m_homogen_dealiasing, 
                                                m_graph, 
                                                m_session->GetVariable(i));
                                    m_base[i]->SetWaveSpace(true);
                                }
                            }
                            else
                            {
                                const LibUtilities::PointsKey PkeyZ(m_npointsZ,
                                        LibUtilities::eFourierEvenlySpaced);
                                const LibUtilities::BasisKey  BkeyZ(
                                        LibUtilities::eFourier,m_npointsZ,
                                        PkeyZ);

                                for (i = 0 ; i < m_base.num_elements(); i++)
                                {
                                    m_base[i] = MemoryManager<MultiRegions
                                        ::ContField3DHomogeneous1D>
                                            ::AllocateSharedPtr(
                                                m_session, BkeyZ, m_LhomZ, 
                                                m_useFFT, m_homogen_dealiasing, 
                                                m_graph, 
                                                m_session->GetVariable(i));
                                    m_base[i]->SetWaveSpace(false);
                                }
                            }
                        }
                        else
                        {
                            i = 0;
                            MultiRegions::ContField2DSharedPtr firstbase =
                                MemoryManager<MultiRegions::ContField2D>
                                ::AllocateSharedPtr(m_session,mesh,
                                                m_session->GetVariable(i));
                            m_base[0]=firstbase;

                            for (i = 1 ; i < m_base.num_elements(); i++)
                            {
                                m_base[i] = MemoryManager<MultiRegions::
                                    ContField2D>::AllocateSharedPtr(
                                        *firstbase,mesh,
                                        m_session->GetVariable(i));
                            }
                        }
                    }
                    break;
                    case 3:
                    {
                        MultiRegions::ContField3DSharedPtr firstbase =
                            MemoryManager<MultiRegions::ContField3D>
                            ::AllocateSharedPtr(m_session, m_graph,
                                                m_session->GetVariable(0));
                        m_base[0] = firstbase;
                        for (i = 1 ; i < m_base.num_elements(); i++)
                        {
                            m_base[i] = MemoryManager<MultiRegions::ContField3D>
                                ::AllocateSharedPtr(*firstbase, m_graph,
                                                    m_session->GetVariable(0));
                        }
                    }
                    break;
                    default:
                        ASSERTL0(false,"Expansion dimension not recognised");
                        break;
                }
            }
            else
            {
                switch(m_expdim)
                {
                    case 1:
                    {
                        // need to use zero for variable as may be more base 
                        // flows than variables
                        for(i = 0 ; i < m_base.num_elements(); i++)
                        {
                            m_base[i] = MemoryManager<MultiRegions
                                ::DisContField1D>
                                ::AllocateSharedPtr(m_session, m_graph,
                                                    m_session->GetVariable(0));
                        }
                        break;
                    }
                    case 2:
                    {
                        for(i = 0 ; i < m_base.num_elements(); i++)
                        {
                            m_base[i] = MemoryManager<MultiRegions::
                                DisContField2D>::AllocateSharedPtr(
                                    m_session,m_graph,
                                    m_session->GetVariable(0));
                        }
                        break;
                    }
                    case 3:
                        ASSERTL0(false, "3 D not set up");
                    default:
                        ASSERTL0(false, "Expansion dimension not recognised");
                        break;
                }
            }
        }
    
        // Import base flow from file and load in m_base        
        void EquationSystem::ImportFldBase(
            std::string pInfile, 
            SpatialDomains::MeshGraphSharedPtr pGraph)
        {
            std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
            std::vector<std::vector<NekDouble> > FieldData;

            //Get Homogeneous
            m_fld->Import(pInfile,FieldDef,FieldData);

            int nvar = m_session->GetVariables().size();
            if (m_session->DefinesSolverInfo("HOMOGENEOUS"))
            {
                std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
            }
            // copy FieldData into m_fields
            for (int j = 0; j < nvar; ++j)
            {
                for(int i = 0; i < FieldDef.size(); ++i)
                {
                    bool flag = FieldDef[i]->m_fields[j] ==
                                m_session->GetVariable(j);
                    ASSERTL0(flag, (std::string("Order of ") + pInfile
                                    + std::string(" data and that defined in "
                                    "the session differs")).c_str());

                    m_base[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
                                                   FieldDef[i]->m_fields[j],
                                                   m_base[j]->UpdateCoeffs());
                }
            }
        }

        /**
         * 
         */
        void EquationSystem::v_DoSolve()
        {

        }

        /**
         * Virtual function to define if operator in DoSolve is
         * negated with regard to the strong form. This is currently
         * only used in Arnoldi solves. Default is false.
         */
        bool EquationSystem::v_NegatedOp(void)
        {
            return false; 
        }

        /**
         * 
         */
        void EquationSystem::v_TransCoeffToPhys()
        {
        
        }
    
        /**
         *
         */
        void EquationSystem::v_TransPhysToCoeff()
        {
        
        }


        /// Virtual function for generating summary information.
        void EquationSystem::v_GenerateSummary(SummaryList& l)
        {
            SessionSummary(l);
        }
        

        /**
         * Write the field data to file. The file is named according to the session
         * name with the extension .fld appended.
         */
        void EquationSystem::v_Output(void)
        {
            WriteFld(m_sessionName + ".fld");
        }
                
        /**
         * Zero the physical fields.
         */
        void EquationSystem::ZeroPhysFields(void)
        {
            for (int i = 0; i < m_fields.num_elements(); i++)
            {
                Vmath::Zero(m_fields[i]->GetNpoints(),
                            m_fields[i]->UpdatePhys(),1);
            }
        }

        /**
         * FwdTrans the m_fields members
         */
        void EquationSystem::FwdTransFields(void)
        {
            for (int i = 0; i < m_fields.num_elements(); i++)
            {
                m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
                                      m_fields[i]->UpdateCoeffs());
                m_fields[i]->SetPhysState(false);
            }
        }

        /**
         * Computes the weak Green form of advection terms (without boundary
         * integral), i.e. \f$ (\nabla \phi \cdot F) \f$ where for example
         * \f$ F=uV \f$.
         * @param   F           Fields.
         * @param   outarray    Storage for result.
         *
         * \note Assuming all fields are of the same expansion and order so that 
         * we can use the parameters of m_fields[0].
         */
        void EquationSystem::WeakAdvectionGreensDivergenceForm(
            const Array<OneD, Array<OneD, NekDouble> > &F,
            Array<OneD, NekDouble> &outarray)
        {
            m_fields[0]->IProductWRTDerivBase(F,outarray);
        }

        /**
         * Calculate Inner product of the divergence advection form
         * \f$(\phi, \nabla \cdot F)\f$, where for example \f$ F = uV \f$.
         * @param   F           Fields.
         * @param   outarray    Storage for result.
         */
        void EquationSystem::WeakAdvectionDivergenceForm(
            const Array<OneD, Array<OneD, NekDouble> > &F,
            Array<OneD, NekDouble> &outarray)
        {
            // Use dimension of Velocity vector to dictate dimension of operation
            int ndim       = F.num_elements();
            int nPointsTot = m_fields[0]->GetNpoints();
            Array<OneD, NekDouble> tmp(nPointsTot);
            Array<OneD, NekDouble> div(nPointsTot, 0.0);

            // Evaluate the divergence
            for (int i = 0; i < ndim; ++i)
            {
                //m_fields[0]->PhysDeriv(i,F[i],tmp);
                m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[i],F[i],tmp);
                Vmath::Vadd(nPointsTot, tmp, 1, div, 1, div, 1);
            }

            m_fields[0]->IProductWRTBase(div, outarray);
        }

        /**
         * Calculate Inner product of the divergence advection form
         * \f$ (\phi, V\cdot \nabla u) \f$
         * @param   V           Fields.
         * @param   u           Fields.
         * @param   outarray    Storage for result.
         */
        void EquationSystem::WeakAdvectionNonConservativeForm(
            const Array<OneD, Array<OneD, NekDouble> > &V,
            const Array<OneD, const NekDouble> &u,
            Array<OneD, NekDouble> &outarray,
            bool UseContCoeffs)
        {
            // use dimension of Velocity vector to dictate dimension of operation
            int ndim       = V.num_elements();

            int nPointsTot = m_fields[0]->GetNpoints();
            Array<OneD, NekDouble> tmp(nPointsTot);
            Array<OneD, NekDouble> wk(ndim * nPointsTot, 0.0);

            AdvectionNonConservativeForm(V, u, tmp, wk);
        
            if (UseContCoeffs)
            {
                m_fields[0]->IProductWRTBase(tmp, outarray,
                                             MultiRegions::eGlobal);
            }
            else
            {
                m_fields[0]->IProductWRTBase_IterPerExp(tmp, outarray);
            }
        }

        /**
         * Calculate the inner product \f$ V\cdot \nabla u \f$
         * @param   V           Fields.
         * @param   u           Fields.
         * @param   outarray    Storage for result.
         * @param   wk          Workspace.
         */
        void EquationSystem::AdvectionNonConservativeForm(
            const Array<OneD, Array<OneD, NekDouble> > &V,
            const Array<OneD, const NekDouble> &u,
            Array<OneD, NekDouble> &outarray,
            Array<OneD, NekDouble> &wk)
        {
            // Use dimension of Velocity vector to dictate dimension of operation
            int ndim       = V.num_elements();
            //int ndim = m_expdim;

            // ToDo: here we should add a check that V has right dimension

            int nPointsTot = m_fields[0]->GetNpoints();
            Array<OneD, NekDouble> grad0,grad1,grad2;

            // Check to see if wk space is defined
            if (wk.num_elements())
            {
                grad0 = wk;
            }
            else
            {
                grad0 = Array<OneD, NekDouble> (nPointsTot);
            }

            // Evaluate V\cdot Grad(u)
            switch(ndim)
            {
                case 1:
                    m_fields[0]->PhysDeriv(u,grad0);
                    Vmath::Vmul(nPointsTot, grad0, 1, V[0], 1, outarray,1);
                    break;
                case 2:
                    grad1 = Array<OneD, NekDouble> (nPointsTot);
                    m_fields[0]->PhysDeriv(u, grad0, grad1);
                    Vmath::Vmul (nPointsTot, grad0, 1, V[0], 1, outarray, 1);
                    Vmath::Vvtvp(nPointsTot, grad1, 1, V[1], 1,
                                 outarray, 1, outarray, 1);
                    break;
                case 3:
                    grad1 = Array<OneD, NekDouble> (nPointsTot);
                    grad2 = Array<OneD, NekDouble> (nPointsTot);
                    m_fields[0]->PhysDeriv(u,grad0,grad1,grad2);
                    Vmath::Vmul (nPointsTot, grad0, 1, V[0], 1, outarray, 1);
                    Vmath::Vvtvp(nPointsTot, grad1, 1, V[1], 1,
                                 outarray, 1, outarray, 1);
                    Vmath::Vvtvp(nPointsTot, grad2, 1, V[2], 1,
                                 outarray, 1, outarray, 1);
                    break;
                default:
                    ASSERTL0(false,"dimension unknown");
            }
        }

        /**
         * @brief Calculate weak DG advection in the form \f$ \langle\phi,
         * \hat{F}\cdot n\rangle - (\nabla \phi \cdot F) \f$
         * 
         * @param   InField                         Fields.
         * @param   OutField                        Storage for result.
         * @param   NumericalFluxIncludesNormal     Default: true.
         * @param   InFieldIsPhysSpace              Default: false.
         * @param   nvariables                      Number of fields.
         */
        void EquationSystem::WeakDGAdvection(
            const Array<OneD, Array<OneD, NekDouble> >& InField,
                  Array<OneD, Array<OneD, NekDouble> >& OutField,
            bool NumericalFluxIncludesNormal,
            bool InFieldIsInPhysSpace,
            int nvariables)
        {
            int i;
            int nVelDim         = m_expdim;
            int nPointsTot      = GetNpoints();
            int ncoeffs         = GetNcoeffs();
            int nTracePointsTot = GetTraceNpoints();

            if (!nvariables)
            {
                nvariables      = m_fields.num_elements();
            }

            Array<OneD, Array<OneD, NekDouble> > fluxvector(nVelDim);
            Array<OneD, Array<OneD, NekDouble> > physfield (nvariables);

            for(i = 0; i < nVelDim; ++i)
            {
                fluxvector[i]    = Array<OneD, NekDouble>(nPointsTot);
            }

            // Get the variables in physical space
            // already in physical space
            if (InFieldIsInPhysSpace == true)
            {
                for (i = 0; i < nvariables; ++i)
                {
                    physfield[i] = InField[i];
                }
            }
            // otherwise do a backward transformation
            else
            {
                for(i = 0; i < nvariables; ++i)
                {
                    // Could make this point to m_fields[i]->UpdatePhys();
                    physfield[i] = Array<OneD, NekDouble>(nPointsTot);
                    m_fields[i]->BwdTrans(InField[i],physfield[i]);
                }
            }

            // Get the advection part (without numerical flux)
            for (i = 0; i < nvariables; ++i)
            {
                // Get the ith component of the  flux vector in (physical space)
                GetFluxVector(i, physfield, fluxvector);

                // Calculate the i^th value of (\grad_i \phi, F)
                WeakAdvectionGreensDivergenceForm(fluxvector,OutField[i]);
            }

            // Get the numerical flux and add to the modal coeffs
            // if the NumericalFluxs function already includes the
            // normal in the output
            if (NumericalFluxIncludesNormal == true)
            {
                Array<OneD, Array<OneD, NekDouble> > numflux   (nvariables);

                for (i = 0; i < nvariables; ++i)
                {
                    numflux[i]   = Array<OneD, NekDouble>(nTracePointsTot);
                }

                // Evaluate numerical flux in physical space which may in
                // general couple all component of vectors
                NumericalFlux(physfield, numflux);

                // Evaulate  <\phi, \hat{F}\cdot n> - OutField[i]
                for (i = 0; i < nvariables; ++i)
                {
                    Vmath::Neg(ncoeffs,OutField[i],1);
                    m_fields[i]->AddTraceIntegral(numflux[i],OutField[i]);
                    m_fields[i]->SetPhysState(false);
                }
            }
            // if the NumericalFlux function does not include the
            // normal in the output
            else
            {
                Array<OneD, Array<OneD, NekDouble> > numfluxX   (nvariables);
                Array<OneD, Array<OneD, NekDouble> > numfluxY   (nvariables);

                for (i = 0; i < nvariables; ++i)
                {
                    numfluxX[i]   = Array<OneD, NekDouble>(nTracePointsTot);
                    numfluxY[i]   = Array<OneD, NekDouble>(nTracePointsTot);
                }

                // Evaluate numerical flux in physical space which may in
                // general couple all component of vectors
                NumericalFlux(physfield, numfluxX, numfluxY);

                // Evaulate  <\phi, \hat{F}\cdot n> - OutField[i]
                for(i = 0; i < nvariables; ++i)
                {
                    Vmath::Neg(ncoeffs,OutField[i],1);
                    m_fields[i]->AddTraceIntegral(numfluxX[i], numfluxY[i],
                                                  OutField[i]);
                    m_fields[i]->SetPhysState(false);
                }
            }
        }
        
        /**
         * Calculate weak DG Diffusion in the LDG form
         * \f$ \langle\psi, \hat{u}\cdot n\rangle
         * - \langle\nabla\psi \cdot u\rangle
         *  \langle\phi, \hat{q}\cdot n\rangle - (\nabla \phi \cdot q) \rangle \f$
         */
        void EquationSystem::WeakDGDiffusion(
            const Array<OneD, Array<OneD, NekDouble> >& InField,
            Array<OneD, Array<OneD, NekDouble> >& OutField,
            bool NumericalFluxIncludesNormal,
            bool InFieldIsInPhysSpace)
        {
            int i, j, k;
            int nPointsTot      = GetNpoints();
            int ncoeffs         = GetNcoeffs();
            int nTracePointsTot = GetTraceNpoints();
            int nvariables      = m_fields.num_elements();
            int nqvar           = 2;

            Array<OneD, NekDouble>  qcoeffs (ncoeffs);
            Array<OneD, NekDouble>  temp (ncoeffs);

            Array<OneD, Array<OneD, NekDouble> > fluxvector (m_spacedim);
            Array<OneD, Array<OneD, NekDouble> > ufield (nvariables);

            Array<OneD, Array<OneD, Array<OneD, NekDouble> > >  flux   (nqvar);
            Array<OneD, Array<OneD, Array<OneD, NekDouble> > >  qfield  (nqvar);

            for (j = 0; j < nqvar; ++j)
            {
                qfield[j]   = Array<OneD, Array<OneD, NekDouble> >(nqvar);
                flux[j]     = Array<OneD, Array<OneD, NekDouble> >(nqvar);

                for (i = 0; i< nvariables; ++i)
                {
                    ufield[i] = Array<OneD, NekDouble>(nPointsTot, 0.0);
                    qfield[j][i]  = Array<OneD, NekDouble>(nPointsTot, 0.0);
                    flux[j][i] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
                }
            }

            for (k = 0; k < m_spacedim; ++k)
            {
                fluxvector[k] = Array<OneD, NekDouble>(nPointsTot, 0.0);
            }

            // Get the variables in physical space already in physical space
            if (InFieldIsInPhysSpace == true)
            {
                for (i = 0; i < nvariables; ++i)
                {
                    ufield[i] = InField[i];
                }
            }
            // Otherwise do a backward transformation
            else
            {
                for (i = 0; i < nvariables; ++i)
                {
                    // Could make this point to m_fields[i]->UpdatePhys();
                    ufield[i] = Array<OneD, NekDouble>(nPointsTot);
                    m_fields[i]->BwdTrans(InField[i],ufield[i]);
                }
            }

            // ##########################################################
            // Compute q_{\eta} and q_{\xi} from su
            // Obtain Numerical Fluxes
            // ##########################################################
            NumFluxforScalar(ufield, flux);

            for (j = 0; j < nqvar; ++j)
            {
                for (i = 0; i < nvariables; ++i)
                {
                    // Get the ith component of the  flux vector in
                    // (physical space) fluxvector = m_tanbasis * u,
                    // where m_tanbasis = 2 by m_spacedim by nPointsTot
                    if (m_tanbasis.num_elements())
                    {
                        for (k = 0; k < m_spacedim; ++k)
                        {
                            Vmath::Vmul(nPointsTot, m_tanbasis[j][k], 1,
                                        ufield[i], 1, fluxvector[k], 1);
                        }
                    }
                    else
                    {
                        GetFluxVector(i, j, ufield, fluxvector);
                    }

                    // Calculate the i^th value of (\grad_i \phi, F)
                    WeakAdvectionGreensDivergenceForm(fluxvector, qcoeffs);

                    Vmath::Neg(ncoeffs,qcoeffs,1);
                    m_fields[i]->AddTraceIntegral(flux[j][i], qcoeffs);
                    m_fields[i]->SetPhysState(false);

                    // Add weighted mass matrix = M ( \nabla \cdot Tanbasis )
//                if(m_gradtan.num_elements())
//                {
//                    MultiRegions::GlobalMatrixKey key(StdRegions::eMass,
//                                                        m_gradtan[j]);
//                    m_fields[i]->MultiRegions::ExpList::GeneralMatrixOp(key,
//                                                        InField[i], temp);
//                    Vmath::Svtvp(ncoeffs, -1.0, temp, 1, qcoeffs, 1,
//                                                        qcoeffs, 1);
//                }

                    //Multiply by the inverse of mass matrix
                    m_fields[i]->MultiplyByElmtInvMass(qcoeffs, qcoeffs);

                    // Back to physical space
                    m_fields[i]->BwdTrans(qcoeffs, qfield[j][i]);
                }
            }


            // ##########################################################
            //   Compute u from q_{\eta} and q_{\xi}
            // ##########################################################

            // Obtain Numerical Fluxes
            NumFluxforVector(ufield, qfield, flux[0]);

            for (i = 0; i < nvariables; ++i)
            {
                // L = L(tan_eta) q_eta + L(tan_xi) q_xi
                OutField[i] = Array<OneD, NekDouble>(ncoeffs, 0.0);
                temp = Array<OneD, NekDouble>(ncoeffs, 0.0);

                if (m_tanbasis.num_elements())
                {
                    for (j = 0; j < nqvar; ++j)
                    {
                        for (k = 0; k < m_spacedim; ++k)
                        {
                            Vmath::Vmul(nPointsTot, m_tanbasis[j][k], 1,
                                        qfield[j][i], 1, fluxvector[k], 1);
                        }

                        WeakAdvectionGreensDivergenceForm(fluxvector, temp);
                        Vmath::Vadd(ncoeffs, temp, 1, OutField[i], 1,
                                    OutField[i], 1);
                    }
                }
                else
                {
                    for (k = 0; k < m_spacedim; ++k)
                    {
                        Vmath::Vcopy(nPointsTot, qfield[k][i], 1,
                                     fluxvector[k], 1);
                    }

                    WeakAdvectionGreensDivergenceForm(fluxvector, OutField[i]);
                }

                // Evaulate  <\phi, \hat{F}\cdot n> - OutField[i]
                Vmath::Neg(ncoeffs,OutField[i],1);
                m_fields[i]->AddTraceIntegral(flux[0][i], OutField[i]);
                m_fields[i]->SetPhysState(false);
            }
        }

        /**
         * Write the n-th checkpoint file.
         * @param   n   The index of the checkpoint file.
         */
        void EquationSystem::Checkpoint_Output(const int n)
        {
            std::string outname =  m_sessionName +  "_" + 
                boost::lexical_cast<std::string>(n);

            WriteFld(outname + ".chk");
        }

        /**
         * Write the n-th checkpoint file.
         * @param   n   The index of the checkpoint file.
         */
        void EquationSystem::Checkpoint_Output(
            const int n, 
            MultiRegions::ExpListSharedPtr &field, 
            std::vector<Array<OneD, NekDouble> > &fieldcoeffs, 
            std::vector<std::string> &variables)
        {
            char chkout[16] = "";
            sprintf(chkout, "%d", n);
            std::string outname = m_sessionName + "_" + chkout + ".chk";
            WriteFld(outname, field, fieldcoeffs, variables);
        }
        
        /**
         * Write the n-th base flow into a .chk file
         * @param   n   The index of the base flow file.
         */
        void EquationSystem::Checkpoint_BaseFlow(const int n)
        {
            std::string outname =  m_sessionName +  "_BaseFlow_" + 
                boost::lexical_cast<std::string>(n);

            WriteFld(outname + ".chk");
        }

        /**
         * Writes the field data to a file with the given filename.
         * @param   outname     Filename to write to.
         */
        void EquationSystem::WriteFld(const std::string &outname)
        {
            std::vector<Array<OneD, NekDouble> > fieldcoeffs(
                m_fields.num_elements());
            std::vector<std::string> variables(m_fields.num_elements());

            for (int i = 0; i < m_fields.num_elements(); ++i)
            {
                if (m_fields[i]->GetNcoeffs() == m_fields[0]->GetNcoeffs())
                {
                    fieldcoeffs[i] = m_fields[i]->UpdateCoeffs();
                }
                else
                {
                    fieldcoeffs[i] = Array<OneD,NekDouble>(m_fields[0]->
                                                           GetNcoeffs());
                    m_fields[0]->ExtractCoeffsToCoeffs(m_fields[i],
                                                       m_fields[i]->GetCoeffs(),
                                                       fieldcoeffs[i]);
                }
                variables[i] = m_boundaryConditions->GetVariable(i);
            }

            v_ExtraFldOutput(fieldcoeffs, variables);

            WriteFld(outname, m_fields[0], fieldcoeffs, variables);
        }



        /**
         * Writes the field data to a file with the given filename.
         * @param   outname         Filename to write to.
         * @param   field           ExpList on which data is based.
         * @param   fieldcoeffs     An array of array of expansion coefficients.
         * @param   variables       An array of variable names.
         */
        void EquationSystem::WriteFld(
            const std::string                    &outname,
            MultiRegions::ExpListSharedPtr       &field,
            std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
            std::vector<std::string>             &variables)
        {
            std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
                = field->GetFieldDefinitions();
            std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());

            // Copy Data into FieldData and set variable
            for(int j = 0; j < fieldcoeffs.size(); ++j)
            {
                for(int i = 0; i < FieldDef.size(); ++i)
                {
                    // Could do a search here to find correct variable
                    FieldDef[i]->m_fields.push_back(variables[j]);
                    field->AppendFieldData(FieldDef[i], FieldData[i],
                                           fieldcoeffs[j]);
                }            
            }

            // Update time in field info if required
            if(m_fieldMetaDataMap.find("Time") != m_fieldMetaDataMap.end())
            {
                m_fieldMetaDataMap["Time"] = boost::lexical_cast<std::string>(m_time);
            }
            
            // If necessary, add mapping information to metadata
            //      and output mapping coordinates
            Array<OneD, MultiRegions::ExpListSharedPtr> fields(1);
            fields[0] = field;           
            GlobalMapping::MappingSharedPtr mapping = 
                    GlobalMapping::Mapping::Load(m_session, fields);
            LibUtilities::FieldMetaDataMap fieldMetaDataMap(m_fieldMetaDataMap);
            mapping->Output( fieldMetaDataMap, outname);

            m_fld->Write(outname, FieldDef, FieldData, fieldMetaDataMap);
        }


        /**
         * Import field from infile and load into \a m_fields. This routine will
         * also perform a \a BwdTrans to ensure data is in both the physical and
         * coefficient storage.
         * @param   infile  Filename to read.
         */
        void EquationSystem::ImportFld(
            const std::string &infile, 
            Array<OneD, MultiRegions::ExpListSharedPtr> &pFields)
        {
            std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
            std::vector<std::vector<NekDouble> > FieldData;

            m_fld->Import(infile,FieldDef,FieldData);

            // Copy FieldData into m_fields
            for(int j = 0; j < pFields.num_elements(); ++j)
            {
                Vmath::Zero(pFields[j]->GetNcoeffs(),
                            pFields[j]->UpdateCoeffs(),1);
                
                for(int i = 0; i < FieldDef.size(); ++i)
                {
                    ASSERTL1(FieldDef[i]->m_fields[j] ==
                             m_session->GetVariable(j),
                             std::string("Order of ") + infile
                             + std::string(" data and that defined in "
                                           "m_boundaryconditions differs"));

                    pFields[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
                                                    FieldDef[i]->m_fields[j],
                                                    pFields[j]->UpdateCoeffs());
                }
                pFields[j]->BwdTrans(pFields[j]->GetCoeffs(),
                                     pFields[j]->UpdatePhys());
            }
        }



        /**
         * Import field from infile and load into \a m_fields. This routine will
         * also perform a \a BwdTrans to ensure data is in both the physical and
         * coefficient storage.
         * @param   infile  Filename to read.
         * If optionan \a ndomains is specified it assumes we loop over nodmains for each nvariables. 
         */
        void EquationSystem::ImportFldToMultiDomains(
                                                     const std::string &infile, 
                                                     Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
                                                     const int ndomains)
        {
            std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
            std::vector<std::vector<NekDouble> > FieldData;
            
            LibUtilities::Import(infile,FieldDef,FieldData);
            
            int nvariables = GetNvariables();

            ASSERTL0(ndomains*nvariables == pFields.num_elements(),"Number of fields does not match the number of variables and domains");
            
            // Copy FieldData into m_fields
            for(int j = 0; j < ndomains; ++j)
            {
                for(int i = 0; i < nvariables; ++i)
                {
                    Vmath::Zero(pFields[j*nvariables+i]->GetNcoeffs(),pFields[j*nvariables+i]->UpdateCoeffs(),1);
                    
                    for(int n = 0; n < FieldDef.size(); ++n)
                    {
                        ASSERTL1(FieldDef[n]->m_fields[i] == m_session->GetVariable(i),
                                 std::string("Order of ") + infile
                                 + std::string(" data and that defined in "
                                               "m_boundaryconditions differs"));
                        
                        pFields[j*nvariables+i]->ExtractDataToCoeffs(FieldDef[n], FieldData[n],
                                                                     FieldDef[n]->m_fields[i],
                                                                     pFields[j*nvariables+i]->UpdateCoeffs());
                    }
                    pFields[j*nvariables+i]->BwdTrans(pFields[j*nvariables+i]->GetCoeffs(),
                                                      pFields[j*nvariables+i]->UpdatePhys());
                }
            }
        }

        /**
         * Import field from infile and load into \a pField. This routine will
         * also perform a \a BwdTrans to ensure data is in both the physical and
         * coefficient storage.
         */
        void EquationSystem::ImportFld(
            const std::string &infile, 
            MultiRegions::ExpListSharedPtr &pField, 
            std::string &pFieldName)
        {
            std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
            std::vector<std::vector<NekDouble> > FieldData;

            m_fld->Import(infile,FieldDef,FieldData);
            int idx = -1;

            Vmath::Zero(pField->GetNcoeffs(),pField->UpdateCoeffs(),1);

            for(int i = 0; i < FieldDef.size(); ++i)
            {
                // find the index of the required field in the file.
                for(int j = 0; j < FieldData.size(); ++j)
                {
                    if (FieldDef[i]->m_fields[j] == pFieldName)
                    {
                        idx = j;
                    }
                }
                ASSERTL1(idx >= 0, "Field " + pFieldName + " not found.");

                pField->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
                                            FieldDef[i]->m_fields[idx],
                                            pField->UpdateCoeffs());
            }
            pField->BwdTrans(pField->GetCoeffs(), pField->UpdatePhys());
        }

        /**
         * Import field from infile and load into the array \a coeffs. 
         *
         * @param infile Filename to read.
         * @param fieldStr an array of string identifying fields to be imported
         * @param coeffs and array of array of coefficients to store imported data
         */
        void EquationSystem::ImportFld(
            const std::string &infile, 
            std::vector< std::string> &fieldStr, 
            Array<OneD, Array<OneD, NekDouble> > &coeffs)
        {

            ASSERTL0(fieldStr.size() <= coeffs.num_elements(),
                     "length of fieldstr should be the same as pFields");
        
            std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
            std::vector<std::vector<NekDouble> > FieldData;
        
            m_fld->Import(infile,FieldDef,FieldData);

            // Copy FieldData into m_fields
            for(int j = 0; j < fieldStr.size(); ++j)
            {
                Vmath::Zero(coeffs[j].num_elements(),coeffs[j],1);
                for(int i = 0; i < FieldDef.size(); ++i)
                {
                    m_fields[0]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
                                                     fieldStr[j], coeffs[j]);
                }
            }
        }

        /**
         * Write out a summary of the session data.
         * @param   out         Output stream to write data to.
         */
        void EquationSystem::SessionSummary(SummaryList& s)
        {
            AddSummaryItem(s, "EquationType", m_session->GetSolverInfo("EQTYPE"));
            AddSummaryItem(s, "Session Name", m_sessionName);
            AddSummaryItem(s, "Spatial Dim.", m_spacedim);
            AddSummaryItem(s, "Max SEM Exp. Order", m_fields[0]->EvalBasisNumModesMax());
            if(m_HomogeneousType == eHomogeneous1D)
            {
                AddSummaryItem(s, "Quasi-3D", "Homogeneous in z-direction");
                AddSummaryItem(s, "Expansion Dim.", m_expdim + 1);
                AddSummaryItem(s, "Num. Hom. Modes (z)", m_npointsZ);
                AddSummaryItem(s, "Hom. length (LZ)", "m_LhomZ");
                AddSummaryItem(s, "FFT Type", m_useFFT ? "FFTW" : "MVM");
                if (m_halfMode)
                {
                    AddSummaryItem(s, "ModeType", "Half Mode");
                }
                else if (m_singleMode)
                {
                    AddSummaryItem(s, "ModeType", "Single Mode");
                }
                else if (m_multipleModes)
                {
                    AddSummaryItem(s, "ModeType", "Multiple Modes");
                    AddSummaryItem(s, "Selected Mode",
                                      boost::lexical_cast<string>(m_NumMode));
                }
            }
            else if(m_HomogeneousType == eHomogeneous2D)
            {
                AddSummaryItem(s, "Quasi-3D", "Homogeneous in yz-plane");
                AddSummaryItem(s, "Expansion Dim.", m_expdim + 2);
                AddSummaryItem(s, "Num. Hom. Modes (y)", m_npointsY);
                AddSummaryItem(s, "Num. Hom. Modes (z)", m_npointsZ);
                AddSummaryItem(s, "Hom. length (LY)", "m_LhomY");
                AddSummaryItem(s, "Hom. length (LZ)", "m_LhomZ");
                AddSummaryItem(s, "FFT Type", m_useFFT ? "FFTW" : "MVM");
            }
            else
            {
                AddSummaryItem(s, "Expansion Dim.", m_expdim);
            }
            
            if (m_session->DefinesSolverInfo("UpwindType"))
            {
                AddSummaryItem(s, "Riemann Solver",
                                  m_session->GetSolverInfo("UpwindType"));
            }
            
            if (m_session->DefinesSolverInfo("AdvectionType"))
            {
                std::string AdvectionType;
                AdvectionType = m_session->GetSolverInfo("AdvectionType");
                AddSummaryItem(s, "Advection Type", GetAdvectionFactory().
                               GetClassDescription(AdvectionType));
            }

            if (m_projectionType == MultiRegions::eGalerkin)
            {
                AddSummaryItem(s, "Projection Type", "Continuous Galerkin");
            }
            else if (m_projectionType == MultiRegions::eDiscontinuous)
            {
                AddSummaryItem(s, "Projection Type", "Discontinuous Galerkin");
            }
            else if (m_projectionType == MultiRegions::eMixed_CG_Discontinuous)
            {
                AddSummaryItem(s, "Projection Type",
                                  "Mixed Continuous Galerkin and Discontinuous");
            }
            
            if (m_session->DefinesSolverInfo("DiffusionType"))
            {
                std::string DiffusionType;
                DiffusionType = m_session->GetSolverInfo("DiffusionType");
                AddSummaryItem(s, "Diffusion Type", GetDiffusionFactory().
                               GetClassDescription(DiffusionType));
            }
        }

        /**
         * Performs a case-insensitive string comparison (from web).
         * @param   s1      First string to compare.
         * @param   s2      Second string to compare.
         * @returns         0 if the strings match.
         */
        int EquationSystem::NoCaseStringCompare(
            const string & s1,
            const string& s2)
        {
            //if (s1.size() < s2.size()) return -1;
            //if (s1.size() > s2.size()) return 1;

            string::const_iterator it1=s1.begin();
            string::const_iterator it2=s2.begin();

            // Stop when either string's end has been reached
            while ( (it1!=s1.end()) && (it2!=s2.end()) )
            {
                if(::toupper(*it1) != ::toupper(*it2)) //letters differ?
                {
                    // Return -1 to indicate smaller than, 1 otherwise
                    return (::toupper(*it1)  < ::toupper(*it2)) ? -1 : 1;
                }

                // Proceed to the next character in each string
                ++it1;
                ++it2;
            }

            size_t size1=s1.size();
            size_t size2=s2.size();// cache lengths

            // Return -1, 0 or 1 according to strings' lengths
            if (size1==size2)
            {
                return 0;
            }

            return (size1 < size2) ? -1 : 1;
        }

        Array<OneD, bool> EquationSystem::v_GetSystemSingularChecks()
        {
            return Array<OneD, bool>(m_session->GetVariables().size(), false);
        }

        void EquationSystem::v_GetFluxVector(
            const int i, Array<OneD,
            Array<OneD, NekDouble> > &physfield,
            Array<OneD, Array<OneD, NekDouble> > &flux)
        {
            ASSERTL0(false, "v_GetFluxVector: This function is not valid "
                     "for the Base class");
        }

        void EquationSystem::v_GetFluxVector(
            const int i, const int j,
            Array<OneD, Array<OneD, NekDouble> > &physfield,
            Array<OneD, Array<OneD, NekDouble> > &flux)
        {
            ASSERTL0(false, "v_GetqFluxVector: This function is not valid "
                     "for the Base class");
        }

        void EquationSystem::v_GetFluxVector(
            const int i, Array<OneD,
            Array<OneD, NekDouble> > &physfield,
            Array<OneD, Array<OneD, NekDouble> > &fluxX,
            Array<OneD, Array<OneD, NekDouble> > &fluxY)
        {
            ASSERTL0(false, "v_GetFluxVector: This function is not valid "
                     "for the Base class");
        }

        void EquationSystem::v_NumericalFlux(
            Array<OneD, Array<OneD, NekDouble> > &physfield,
            Array<OneD, Array<OneD, NekDouble> > &numflux)
        {
            ASSERTL0(false, "v_NumericalFlux: This function is not valid "
                     "for the Base class");
        }

        void EquationSystem::v_NumericalFlux(
            Array<OneD, Array<OneD, NekDouble> > &physfield,
            Array<OneD, Array<OneD, NekDouble> > &numfluxX,
            Array<OneD, Array<OneD, NekDouble> > &numfluxY )
        {
            ASSERTL0(false, "v_NumericalFlux: This function is not valid "
                     "for the Base class");
        }

        void EquationSystem::v_NumFluxforScalar(
            const Array<OneD, Array<OneD, NekDouble> >         &ufield,
            Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &uflux)
        {
            ASSERTL0(false, "v_NumFluxforScalar: This function is not valid "
                     "for the Base class");
        }

        void EquationSystem::v_NumFluxforVector(
            const Array<OneD, Array<OneD, NekDouble> >   &ufield,
            Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &qfield,
            Array<OneD, Array<OneD, NekDouble > >              &qflux)
        {
            ASSERTL0(false, "v_NumFluxforVector: This function is not valid "
                     "for the Base class");
        }

        MultiRegions::ExpListSharedPtr EquationSystem::v_GetPressure()
        {
            ASSERTL0(false, "This function is not valid for the Base class");
            MultiRegions::ExpListSharedPtr null;
            return null;
        }

        void EquationSystem::v_ExtraFldOutput(
            std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
            std::vector<std::string>             &variables)
        {
        }

    }
}