MeshGraph.cpp

////////////////////////////////////////////////////////////////////////////////
//
//  File:  $Source: /usr/sci/projects/Nektar/cvs/Nektar++/library/SpatialDomains/MeshGraph.cpp,v $
//
//  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:
//
//
////////////////////////////////////////////////////////////////////////////////
#include "pchSpatialDomains.h"

#include <boost/foreach.hpp>

#include <SpatialDomains/MeshGraph.h>
#include <LibUtilities/BasicUtils/ParseUtils.hpp>

// Use the stl version, primarily for string.
#ifndef TIXML_USE_STL
#define TIXML_USE_STL
#endif

#include <tinyxml/tinyxml.h>
#include <cstring>
#include <sstream>
#include <cmath>

#include <SpatialDomains/MeshGraph1D.h>
#include <SpatialDomains/MeshGraph2D.h>
#include <SpatialDomains/MeshGraph3D.h>

// These are required for the Write(...) and Import(...) functions.
#include <boost/archive/iterators/base64_from_binary.hpp>
#include <boost/archive/iterators/binary_from_base64.hpp>
#include <boost/archive/iterators/transform_width.hpp>
#include <boost/iostreams/copy.hpp>
#include <boost/iostreams/filter/zlib.hpp>
#include <boost/iostreams/filtering_stream.hpp>

namespace Nektar
{
    namespace SpatialDomains
    {

        MeshGraph::MeshGraph():
            m_meshDimension(3),
            m_spaceDimension(3)
        {
        }

        MeshGraph::MeshGraph(unsigned int meshDimension, unsigned int spaceDimension) :
            m_meshDimension(meshDimension),
            m_spaceDimension(spaceDimension)
        {
        }

        MeshGraph::MeshGraph(const LibUtilities::SessionReaderSharedPtr &pSession)
            : m_session(pSession)
        {

        }

        boost::shared_ptr<MeshGraph> MeshGraph::Read(const std::string& infilename, bool pReadExpansions)
        {
            boost::shared_ptr<MeshGraph> returnval;

            MeshGraph mesh;

            mesh.ReadGeometry(infilename);
            int meshDim = mesh.GetMeshDimension();

            switch(meshDim)
            {
            case 1:
                returnval = MemoryManager<MeshGraph1D>::AllocateSharedPtr();
                break;

            case 2:
                returnval = MemoryManager<MeshGraph2D>::AllocateSharedPtr();
                break;

            case 3:
                returnval = MemoryManager<MeshGraph3D>::AllocateSharedPtr();
                break;

            default:
                std::string err = "Invalid mesh dimension: ";
                std::stringstream strstrm;
                strstrm << meshDim;
                err += strstrm.str();
                NEKERROR(ErrorUtil::efatal, err.c_str());
            }

            if (returnval)
            {
                returnval->ReadGeometry(infilename);
                returnval->ReadGeometryInfo(infilename);
                if (pReadExpansions)
                {
                    returnval->ReadExpansions(infilename);
                }
            }
            return returnval;
        }

        void MeshGraph::SetExpansions(std::vector<SpatialDomains::FieldDefinitionsSharedPtr> &fielddef)
        {
            int i,j,k,g,h,cnt,id;
            GeometrySharedPtr geom;
            
            ExpansionMapShPtr expansionMap;

            // Loop over fields and determine unique fields string and
            // declare whole expansion list
            for(i = 0; i < fielddef.size(); ++i)
            {
                for(j = 0; j < fielddef[i]->m_fields.size(); ++j)
                {
                    std::string field = fielddef[i]->m_fields[j];
                    if(m_expansionMapShPtrMap.count(field) == 0)
                    {
                        expansionMap = MemoryManager<ExpansionMap>::AllocateSharedPtr();
                        m_expansionMapShPtrMap[field] = expansionMap;
                        
                        // check to see if DefaultVar also not set and if so assign it to this expansion
                        if(m_expansionMapShPtrMap.count("DefaultVar") == 0)
                        {
                            m_expansionMapShPtrMap["DefaultVar"] = expansionMap;
                        }

                        // loop over all elements and set expansion
                        for(k = 0; k < fielddef.size(); ++k)
                        {
                            for(h = 0; h < fielddef[k]->m_fields.size(); ++h)
                            {
                                if(fielddef[k]->m_fields[h] == field)
                                {
                                    expansionMap = m_expansionMapShPtrMap.find(field)->second;
                                    LibUtilities::BasisKeyVector def;
                                    
                                    for(g = 0; g < fielddef[k]->m_elementIDs.size(); ++g)
                                    {
                                        ExpansionShPtr tmpexp =
                                            MemoryManager<Expansion>::AllocateSharedPtr(geom, def);
                                        (*expansionMap)[fielddef[k]->m_elementIDs[g]] = tmpexp;
                                    }
                                }
                            }
                        }
                    }
                }
            }                    
            
            // loop over all elements find the geometry shared ptr and
            // set up basiskey vector
            for(i = 0; i < fielddef.size(); ++i)
            {
                cnt = 0;
                std::vector<std::string>  fields = fielddef[i]->m_fields;
                std::vector<unsigned int> nmodes = fielddef[i]->m_numModes;
                std::vector<LibUtilities::BasisType> basis = fielddef[i]->m_basis;
                bool pointDef = fielddef[i]->m_pointsDef;
                bool numPointDef = fielddef[i]->m_numPointsDef;

                // Check points and numpoints
                std::vector<unsigned int> npoints = fielddef[i]->m_numPoints;
                std::vector<LibUtilities::PointsType> points = fielddef[i]->m_points;

                bool UniOrder =  fielddef[i]->m_uniOrder;

                int check = 0;
                for(j=0; j< basis.size(); ++j)
                {
                    if( (strcmp(LibUtilities::BasisTypeMap[basis[j]], "Modified_A") == 0) ||
                        (strcmp(LibUtilities::BasisTypeMap[basis[j]], "Modified_B") == 0) ||
                        (strcmp(LibUtilities::BasisTypeMap[basis[j]], "Modified_C") == 0) ||
                        (strcmp(LibUtilities::BasisTypeMap[basis[j]], "Fourier") == 0) )
                    {
                        check++;
                    }
                }

                if(check==basis.size())
                {
                    for(j = 0; j < fielddef[i]->m_elementIDs.size(); ++j)
                    {

                        LibUtilities::BasisKeyVector bkeyvec;
                        id = fielddef[i]->m_elementIDs[j];
                        
                        switch(fielddef[i]->m_shapeType)
                        {
                        case eSegment:
                            {
                                ASSERTL0(m_segGeoms.count(fielddef[i]->m_elementIDs[j]),
                                        "Failed to find geometry with same global id");
                                geom = m_segGeoms[fielddef[i]->m_elementIDs[j]];

                                LibUtilities::PointsKey pkey(nmodes[cnt]+1,LibUtilities::eGaussLobattoLegendre);

                                if(numPointDef&&pointDef)
                                {
                                    const LibUtilities::PointsKey pkey1(npoints[cnt],points[0]);
                                    pkey = pkey1;
                                }
                                else if(!numPointDef&&pointDef)
                                {
                                    const LibUtilities::PointsKey pkey1(nmodes[cnt]+1,points[0]);
                                    pkey = pkey1;
                                }
                                else if(numPointDef&&!pointDef)
                                {
                                    const LibUtilities::PointsKey pkey1(npoints[cnt],LibUtilities::eGaussLobattoLegendre);
                                    pkey = pkey1;
                                }

                                LibUtilities::BasisKey bkey(basis[0],nmodes[cnt],pkey);

                                if(!UniOrder)
                                {
                                    cnt++;
                                }
                                bkeyvec.push_back(bkey);
                            }
                            break;
                        case eTriangle:
                            {
                                ASSERTL0(m_triGeoms.count(fielddef[i]->m_elementIDs[j]),
                                        "Failed to find geometry with same global id");
                                geom = m_triGeoms[fielddef[i]->m_elementIDs[j]];

                                LibUtilities::PointsKey pkey(nmodes[cnt]+1,LibUtilities::eGaussLobattoLegendre);
                                if(numPointDef&&pointDef)
                                {
                                    const LibUtilities::PointsKey pkey2(npoints[cnt],points[0]);
                                    pkey = pkey2;
                                }
                                else if(!numPointDef&&pointDef)
                                {
                                    const LibUtilities::PointsKey pkey2(nmodes[cnt]+1,points[0]);
                                    pkey = pkey2;
                                }
                                else if(numPointDef&&!pointDef)
                                {
                                    const LibUtilities::PointsKey pkey2(npoints[cnt],LibUtilities::eGaussLobattoLegendre);
                                    pkey = pkey2;
                                }
                                LibUtilities::BasisKey bkey(basis[0],nmodes[cnt],pkey);

                                bkeyvec.push_back(bkey);

                                LibUtilities::PointsKey pkey1(nmodes[cnt+1],LibUtilities::eGaussRadauMAlpha1Beta0);
                                if(numPointDef&&pointDef)
                                {
                                    const LibUtilities::PointsKey pkey2(npoints[cnt+1],points[1]);
                                    pkey1 = pkey2;
                                }
                                else if(!numPointDef&&pointDef)
                                {
                                    const LibUtilities::PointsKey pkey2(nmodes[cnt+1],points[1]);
                                    pkey1 = pkey2;
                                }
                                else if(numPointDef&&!pointDef)
                                {
                                    const LibUtilities::PointsKey pkey2(npoints[cnt+1],LibUtilities::eGaussRadauMAlpha1Beta0);
                                    pkey1 = pkey2;
                                }
                                LibUtilities::BasisKey bkey1(basis[1],nmodes[cnt+1],pkey1);
                                bkeyvec.push_back(bkey1);

                                if(!UniOrder)
                                {
                                    cnt += 2;
                                }
                            }
                            break;
                        case eQuadrilateral:
                            {
                                ASSERTL0(m_quadGeoms.count(fielddef[i]->m_elementIDs[j]),
                                        "Failed to find geometry with same global id");
                                geom = m_quadGeoms[fielddef[i]->m_elementIDs[j]];

                                for(int b = 0; b < 2; ++b)
                                {
                                    LibUtilities::PointsKey pkey(nmodes[cnt+b]+1,LibUtilities::eGaussLobattoLegendre);

                                    if(numPointDef&&pointDef)
                                    {
                                        const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
                                        pkey = pkey2;
                                    }
                                    else if(!numPointDef&&pointDef)
                                    {
                                        const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
                                        pkey = pkey2;
                                    }
                                    else if(numPointDef&&!pointDef)
                                    {
                                        const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
                                        pkey = pkey2;
                                    }
                                    LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
                                    bkeyvec.push_back(bkey);
                                }

                                if(!UniOrder)
                                {
                                    cnt += 2;
                                }
                            }
                            break;

                    case eTetrahedron:
                        {
                            k = fielddef[i]->m_elementIDs[j];
                            ASSERTL0(m_tetGeoms.find(k) != m_tetGeoms.end(),
                                    "Failed to find geometry with same global id");
                            geom = m_tetGeoms[k];
                            
                            for(int b = 0; b < 3; ++b)
                            {
                                LibUtilities::PointsKey pkey(nmodes[cnt+b],points[b]);

                                if(numPointDef&&pointDef)
                                {
                                    const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
                                    pkey = pkey2;
                                }
                                else if(!numPointDef&&pointDef)
                                {
                                    const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
                                    pkey = pkey2;
                                }
                                else if(numPointDef&&!pointDef)
                                {
                                    const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
                                    pkey = pkey2;
                                }
                                
                                LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);

                                bkeyvec.push_back(bkey);
                            }

                            if(!UniOrder)
                            {
                                cnt += 3;
                            }
                        }
                        break;
                    case ePrism:
                        {
                            k = fielddef[i]->m_elementIDs[j];
                            ASSERTL0(m_prismGeoms.find(k) != m_prismGeoms.end(),
                                    "Failed to find geometry with same global id");
                            geom = m_prismGeoms[k];

                            for(int b = 0; b < 3; ++b)
                            {
                                LibUtilities::PointsKey pkey(nmodes[cnt+b],points[b]);
                                
                                if(numPointDef&&pointDef)
                                {
                                    const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
                                    pkey = pkey2;
                                }
                                else if(!numPointDef&&pointDef)
                                {
                                    const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
                                    pkey = pkey2;
                                }
                                else if(numPointDef&&!pointDef)
                                {
                                    const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
                                    pkey = pkey2;
                                }
                            
                                LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
                                bkeyvec.push_back(bkey);
                            }
                            
                            if(!UniOrder)
                            {
                                cnt += 3;
                            }
                        }
                        break;
                    case eHexahedron:
                        {
                            k = fielddef[i]->m_elementIDs[j];
                            ASSERTL0(m_hexGeoms.find(k) != m_hexGeoms.end(),
                                    "Failed to find geometry with same global id");
                            geom = m_hexGeoms[k];

                            for(int b = 0; b < 3; ++b)
                            {
                                LibUtilities::PointsKey pkey(nmodes[cnt+b],points[b]);

                                if(numPointDef&&pointDef)
                                {
                                    const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
                                    pkey = pkey2;
                                }
                                else if(!numPointDef&&pointDef)
                                {
                                    const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
                                    pkey = pkey2;
                                }
                                else if(numPointDef&&!pointDef)
                                {
                                    const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
                                    pkey = pkey2;
                                }
                                
                                LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
                                bkeyvec.push_back(bkey);
                            }

                            if(!UniOrder)
                            {
                                cnt += 3;
                            }
                        }
                        break;
                        default:
                            ASSERTL0(false,"Need to set up for pyramid and prism 3D Expansions");
                            break;
                        }
                        
                        for(k = 0; k < fields.size(); ++k)
                        {
                            expansionMap = m_expansionMapShPtrMap.find(fields[k])->second;
                            (*expansionMap)[id]->m_geomShPtr = geom;
                            (*expansionMap)[id]->m_basisKeyVector = bkeyvec;
                        }
                    }
                }
                else
                {
                    ASSERTL0(false,"Need to set up for non Modified basis");
                }
            }
        }


        void MeshGraph::SetExpansions(std::vector<SpatialDomains::FieldDefinitionsSharedPtr> &fielddef, std::vector< std::vector<LibUtilities::PointsType> > &pointstype)
        {
            
            int i,j,k,g,h,cnt,id;
            GeometrySharedPtr geom;

            ExpansionMapShPtr expansionMap;

            // Loop over fields and determine unique fields string and
            // declare whole expansion list
            for(i = 0; i < fielddef.size(); ++i)
            {
                for(j = 0; j < fielddef[i]->m_fields.size(); ++j)
                {
                    std::string field = fielddef[i]->m_fields[j];
                    if(m_expansionMapShPtrMap.count(field) == 0)
                    {
                        expansionMap = MemoryManager<ExpansionMap>::AllocateSharedPtr();
                        m_expansionMapShPtrMap[field] = expansionMap;
                        
                        // check to see if DefaultVar also not set and if so assign it to this expansion
                        if(m_expansionMapShPtrMap.count("DefaultVar") == 0)
                        {
                            m_expansionMapShPtrMap["DefaultVar"] = expansionMap;
                        }

                        // loop over all elements and set expansion
                        for(k = 0; k < fielddef.size(); ++k)
                        {
                            for(h = 0; h < fielddef[k]->m_fields.size(); ++h)
                            {
                                if(fielddef[k]->m_fields[h] == field)
                                {
                                    expansionMap = m_expansionMapShPtrMap.find(field)->second;
                                    LibUtilities::BasisKeyVector def;
                                    
                                    for(g = 0; g < fielddef[k]->m_elementIDs.size(); ++g)
                                    {
                                        ExpansionShPtr tmpexp =
                                            MemoryManager<Expansion>::AllocateSharedPtr(geom, def);
                                        (*expansionMap)[fielddef[k]->m_elementIDs[g]] = tmpexp;
                                    }
                                }
                            }
                        }
                    }
                }
            }                    


            // loop over all elements find the geometry shared ptr and
            // set up basiskey vector
            for(i = 0; i < fielddef.size(); ++i)
            {
                cnt = 0;
                std::vector<std::string>  fields = fielddef[i]->m_fields;
                std::vector<unsigned int> nmodes = fielddef[i]->m_numModes;
                std::vector<LibUtilities::BasisType> basis = fielddef[i]->m_basis;
                bool UniOrder =  fielddef[i]->m_uniOrder;

                for(j = 0; j < fielddef[i]->m_elementIDs.size(); ++j)
                {

                    LibUtilities::BasisKeyVector bkeyvec;
                    id = fielddef[i]->m_elementIDs[j];

                    switch(fielddef[i]->m_shapeType)
                    {
                    case eSegment:
                        {
                            k = fielddef[i]->m_elementIDs[j];
                            ASSERTL0(m_segGeoms.find(k) != m_segGeoms.end(),
                                     "Failed to find geometry with same global id.");
                            geom = m_segGeoms[k];
//                            for(k = 0; k < m_segGeoms.size();++k)
//                            {
//                                if(m_segGeoms[k]->GetGlobalID() == fielddef[i]->m_elementIDs[j])
//                                {
//                                    geom = m_segGeoms[k];
//                                    break;
//                                }
//                            }
//                            ASSERTL0(k != m_segGeoms.size(),"Failed to find geometry with same global id");

                            const LibUtilities::PointsKey pkey(nmodes[cnt],pointstype[i][0]);
                            LibUtilities::BasisKey bkey(basis[0],nmodes[cnt],pkey);
                            if(!UniOrder)
                            {
                                cnt++;
                            }
                            bkeyvec.push_back(bkey);
                        }
                        break;
                    case eTriangle:
                        {
                            k = fielddef[i]->m_elementIDs[j];
                            ASSERTL0(m_triGeoms.find(k) != m_triGeoms.end(),
                                     "Failed to find geometry with same global id.");
                            geom = m_triGeoms[k];
//                            for(k = 0; k < m_triGeoms.size();++k)
//                            {
//                                if(m_triGeoms[k]->GetGlobalID() == fielddef[i]->m_elementIDs[j])
//                                {
//                                    geom = m_triGeoms[k];
//                                    break;
//                                }
//                            }
//                            ASSERTL0(k != m_triGeoms.size(),"Failed to find geometry with same global id");
                            for(int b = 0; b < 2; ++b)
                            {
                                const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
                                LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
                                bkeyvec.push_back(bkey);
                            }

                            if(!UniOrder)
                            {
                                cnt += 2;
                            }
                        }
                        break;
                    case eQuadrilateral:
                        {
                            k = fielddef[i]->m_elementIDs[j];
                            ASSERTL0(m_quadGeoms.find(k) != m_quadGeoms.end(),
                                    "Failed to find geometry with same global id");
                            geom = m_quadGeoms[k];
//                            for(k = 0; k < m_quadGeoms.size();++k)
//                            {
//                                if(m_quadGeoms[k]->GetGlobalID() == fielddef[i]->m_elementIDs[j])
//                                {
//                                    geom = m_quadGeoms[k];
//                                    break;
//                                }
//                            }
//                            ASSERTL0(k != m_quadGeoms.size(),"Failed to find geometry with same global id");

                            for(int b = 0; b < 2; ++b)
                            {
                                const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
                                LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
                                bkeyvec.push_back(bkey);
                            }

                            if(!UniOrder)
                            {
                                cnt += 2;
                            }
                        }
                        break;
                    case eTetrahedron:
                        {
                            k = fielddef[i]->m_elementIDs[j];
                            ASSERTL0(m_tetGeoms.find(k) != m_tetGeoms.end(),
                                    "Failed to find geometry with same global id");
                            geom = m_tetGeoms[k];

//                            for(k = 0; k < m_tetGeoms.size();++k)
//                            {
//                                if(m_tetGeoms[k]->GetGlobalID() == fielddef[i]->m_elementIDs[j])
//                                {
//                                    geom = m_tetGeoms[k];
//                                    break;
//                                }
//                            }
//                            ASSERTL0(k != m_tetGeoms.size(),"Failed to find geometry with same global id");
                            for(int b = 0; b < 3; ++b)
                            {
                                const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
                                LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
                                bkeyvec.push_back(bkey);
                            }

                            if(!UniOrder)
                            {
                                cnt += 2;
                            }
                        }
                        break;
                    case ePrism:
                        {
                            k = fielddef[i]->m_elementIDs[j];
                            ASSERTL0(m_prismGeoms.find(k) != m_prismGeoms.end(),
                                    "Failed to find geometry with same global id");
                            geom = m_prismGeoms[k];

                            for(int b = 0; b < 3; ++b)
                            {
                                const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
                                LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
                                bkeyvec.push_back(bkey);
                            }

                            if(!UniOrder)
                            {
                                cnt += 2;
                            }
                        }
                        break;
                    case eHexahedron:
                        {
                            k = fielddef[i]->m_elementIDs[j];
                            ASSERTL0(m_hexGeoms.find(k) != m_hexGeoms.end(),
                                    "Failed to find geometry with same global id");
                            geom = m_hexGeoms[k];
//                            for(k = 0; k < m_quadGeoms.size();++k)
//                            {
//                                if(m_hexGeoms[k]->GetGlobalID() == fielddef[i]->m_elementIDs[j])
//                                {
//                                    geom = m_hexGeoms[k];
//                                    break;
//                                }
//                            }
//                            ASSERTL0(k != m_hexGeoms.size(),"Failed to find geometry with same global id");

                            for(int b = 0; b < 3; ++b)
                            {
                                const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
                                LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
                                bkeyvec.push_back(bkey);
                            }

                            if(!UniOrder)
                            {
                                cnt += 2;
                            }
                        }
                        break;
                    default:
                        ASSERTL0(false,"Need to set up for pyramid and prism 3D Expansions");
                        break;
                    }

                        for(k = 0; k < fields.size(); ++k)
                        {
                            expansionMap = m_expansionMapShPtrMap.find(fields[k])->second;
                            (*expansionMap)[id]->m_geomShPtr = geom;
                            (*expansionMap)[id]->m_basisKeyVector = bkeyvec;
                        }
                }
            }

        }


        // \brief Read will read the meshgraph vertices given a filename.
        void MeshGraph::ReadGeometry(const std::string& infilename)
        {
          TiXmlDocument doc(infilename);
          bool loadOkay = doc.LoadFile();

          std::stringstream errstr;
          errstr << "Unable to load file: " << infilename << " (";
          errstr << doc.ErrorDesc() << ", line " << doc.ErrorRow()
                 << ", column " << doc.ErrorCol() << ")";
          ASSERTL0(loadOkay, errstr.str());

          ReadGeometry(doc);
        }

        // \brief Read will read the meshgraph vertices given a TiXmlDocument.
        void MeshGraph::ReadGeometry(TiXmlDocument &doc)
        {
            TiXmlHandle docHandle(&doc);
            TiXmlNode* node = NULL;
            TiXmlElement* mesh = NULL;
            TiXmlElement* master = NULL;    // Master tag within which all data is contained.

            int err;    /// Error value returned by TinyXML.

            master = doc.FirstChildElement("NEKTAR");
            ASSERTL0(master, "Unable to find NEKTAR tag in file.");

            // Find the Mesh tag and same the dim and space attributes
            mesh = master->FirstChildElement("GEOMETRY");

            ASSERTL0(mesh, "Unable to find GEOMETRY tag in file.");
            TiXmlAttribute *attr = mesh->FirstAttribute();

            // Initialize the mesh and space dimensions to 3 dimensions.
            // We want to do this each time we read a file, so it should
            // be done here and not just during class initialization.
            m_meshPartitioned = false;
            m_meshDimension = 3;
            m_spaceDimension = 3;

            while (attr)
            {
                std::string attrName(attr->Name());
                if (attrName == "DIM")
                {
                    err = attr->QueryIntValue(&m_meshDimension);
                    ASSERTL1(err==TIXML_SUCCESS, "Unable to read mesh dimension.");
                }
                else if (attrName == "SPACE")
                {
                    err = attr->QueryIntValue(&m_spaceDimension);
                    ASSERTL1(err==TIXML_SUCCESS, "Unable to read space dimension.");
                }
                else if (attrName == "PARTITION")
                {
                    err = attr->QueryIntValue(&m_partition);
                    ASSERTL1(err==TIXML_SUCCESS, "Unable to read partition.");
                    m_meshPartitioned = true;
                }
                else
                {
                    std::string errstr("Unknown attribute: ");
                    errstr += attrName;
                    ASSERTL1(false, errstr.c_str());
                }

                // Get the next attribute.
                attr = attr->Next();
            }

            ASSERTL1(m_meshDimension<=m_spaceDimension, "Mesh dimension greater than space dimension");

            // Now read the vertices
            TiXmlElement* element = mesh->FirstChildElement("VERTEX");
            ASSERTL0(element, "Unable to find mesh VERTEX tag in file.");

            NekDouble xscale,yscale,zscale;
            
            // check to see if any scaling parameters are in
            // attributes and determine these values
            LibUtilities::ExpressionEvaluator expEvaluator;
            const char *xscal =  element->Attribute("XSCALE");
            if(!xscal)
            {
                xscale = 1.0;
            }
            else
            {
                std::string xscalstr = xscal;
                expEvaluator.DefineFunction("",xscalstr);
                xscale = expEvaluator.Evaluate();
            }

            const char *yscal =  element->Attribute("YSCALE");
            if(!yscal)
            {
                yscale = 1.0;
            }
            else
            {
                std::string yscalstr = yscal;
                expEvaluator.DefineFunction("",yscalstr);
                yscale = expEvaluator.Evaluate();
            }
                
            const char *zscal = element->Attribute("ZSCALE");
            if(!zscal)
            {
                zscale = 1.0;
            }
            else
            {
                std::string zscalstr = zscal;
                expEvaluator.DefineFunction("",zscalstr);
                zscale = expEvaluator.Evaluate();
            }
            
            TiXmlElement *vertex = element->FirstChildElement("V");

            int indx;
            int nextVertexNumber = -1;

            while (vertex)
            {
                nextVertexNumber++;

                TiXmlAttribute *vertexAttr = vertex->FirstAttribute();
                std::string attrName(vertexAttr->Name());

                ASSERTL0(attrName == "ID", (std::string("Unknown attribute name: ") + attrName).c_str());

                err = vertexAttr->QueryIntValue(&indx);
                ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
//                ASSERTL0(indx == nextVertexNumber, "Vertex IDs must begin with zero and be sequential.");

                // Now read body of vertex
                std::string vertexBodyStr;

                TiXmlNode *vertexBody = vertex->FirstChild();

                while (vertexBody)
                {
                    // Accumulate all non-comment body data.
                    if (vertexBody->Type() == TiXmlNode::TEXT)
                    {
                        vertexBodyStr += vertexBody->ToText()->Value();
                        vertexBodyStr += " ";
                    }

                    vertexBody = vertexBody->NextSibling();
                }

                ASSERTL0(!vertexBodyStr.empty(), "Vertex definitions must contain vertex data.");

                // Get vertex data from the data string.
                NekDouble xval, yval, zval;
                std::istringstream vertexDataStrm(vertexBodyStr.c_str());

                try
                {
                    while(!vertexDataStrm.fail())
                    {
                        vertexDataStrm >> xval >> yval >> zval;
                        
                        xval *= xscale;
                        yval *= yscale;
                        zval *= zscale;
                        
                        // Need to check it here because we may not be
                        // good after the read indicating that there
                        // was nothing to read.
                        if (!vertexDataStrm.fail())
                        {
                            VertexComponentSharedPtr vert(MemoryManager<VertexComponent>::AllocateSharedPtr(m_spaceDimension, indx, xval, yval, zval));
                            m_vertSet[indx] = vert;
                        }
                    }
                }
                catch(...)
                {
                    ASSERTL0(false, "Unable to read VERTEX data.");
                }

                vertex = vertex->NextSiblingElement("V");
            }
        }


        void MeshGraph::WriteGeometry(const std::string& fileName)
        {
            TiXmlDocument doc(fileName);
            WriteGeometry(doc);
        }

        void MeshGraph::WriteGeometry(TiXmlDocument& doc)
        {

        }

        VertexComponentSharedPtr MeshGraph::AddVertex(NekDouble x, NekDouble y, NekDouble z)
        {
            unsigned int nextId = m_vertSet.rbegin()->first + 1;
            VertexComponentSharedPtr vert(MemoryManager<VertexComponent>::AllocateSharedPtr(m_spaceDimension, nextId, x, y, z));
            m_vertSet[nextId] = vert;
            return vert;
        }

        SegGeomSharedPtr MeshGraph::AddEdge(VertexComponentSharedPtr v0, VertexComponentSharedPtr v1,
            CurveSharedPtr curveDefinition)
        {
            VertexComponentSharedPtr vertices[] = {v0, v1};
            SegGeomSharedPtr edge;
            int edgeId = m_segGeoms.rbegin()->first + 1;

            if( curveDefinition )
            {
                edge = MemoryManager<SegGeom>::AllocateSharedPtr(edgeId, m_spaceDimension, vertices, curveDefinition);
            }
            else
            {
                edge = MemoryManager<SegGeom>::AllocateSharedPtr(edgeId, m_spaceDimension, vertices);
            }
            m_segGeoms[edgeId] = edge;
            return edge;
        }

        TriGeomSharedPtr MeshGraph::AddTriangle(SegGeomSharedPtr edges[], StdRegions::EdgeOrientation orient[])
        {
            int indx = m_triGeoms.rbegin()->first + 1;
            TriGeomSharedPtr trigeom(MemoryManager<TriGeom>::AllocateSharedPtr(indx, edges, orient));
            trigeom->SetGlobalID(indx);

            m_triGeoms[indx] = trigeom;

            return trigeom;
        }

        QuadGeomSharedPtr MeshGraph::AddQuadrilateral(SegGeomSharedPtr edges[], StdRegions::EdgeOrientation orient[])
        {
            int indx = m_quadGeoms.rbegin()->first + 1;
            QuadGeomSharedPtr quadgeom(MemoryManager<QuadGeom>::AllocateSharedPtr(indx, edges, orient));
            quadgeom->SetGlobalID(indx);

            m_quadGeoms[indx] = quadgeom;
            return quadgeom;
        }

        PrismGeomSharedPtr MeshGraph::AddPrism(TriGeomSharedPtr tfaces[PrismGeom::kNtfaces],
            QuadGeomSharedPtr qfaces[PrismGeom::kNqfaces])
        {
            // Setting the orientation is disabled in the reader.  Why?
            StdRegions::FaceOrientation faceorient[PrismGeom::kNtfaces + PrismGeom::kNqfaces];

            unsigned int index = m_prismGeoms.rbegin()->first + 1;
            PrismGeomSharedPtr prismgeom(MemoryManager<PrismGeom>::AllocateSharedPtr(tfaces, qfaces, faceorient));
            prismgeom->SetGlobalID(index);

            m_prismGeoms[index] = prismgeom;
            return prismgeom;
        }

        TetGeomSharedPtr MeshGraph::AddTetrahedron(TriGeomSharedPtr tfaces[TetGeom::kNtfaces])
        {
            unsigned int index = m_tetGeoms.rbegin()->first + 1;
            TetGeomSharedPtr tetgeom(MemoryManager<TetGeom>::AllocateSharedPtr(tfaces));
            tetgeom->SetGlobalID(index);

            m_tetGeoms[index] = tetgeom;
            return tetgeom;
        }

        PyrGeomSharedPtr MeshGraph::AddPyramid(TriGeomSharedPtr tfaces[PyrGeom::kNtfaces],
            QuadGeomSharedPtr qfaces[PyrGeom::kNqfaces])
        {
            unsigned int index = m_pyrGeoms.rbegin()->first + 1;

            StdRegions::FaceOrientation faceorient[PyrGeom::kNtfaces + PyrGeom::kNqfaces];
            PyrGeomSharedPtr pyrgeom(MemoryManager<PyrGeom>::AllocateSharedPtr(tfaces, qfaces, faceorient));
            pyrgeom->SetGlobalID(index);

            m_pyrGeoms[index] = pyrgeom;
            return pyrgeom;
        }

        HexGeomSharedPtr MeshGraph::AddHexahedron(QuadGeomSharedPtr qfaces[HexGeom::kNqfaces])
        {
            unsigned int index = m_hexGeoms.rbegin()->first + 1;
            HexGeomSharedPtr hexgeom(MemoryManager<HexGeom>::AllocateSharedPtr(qfaces));
            hexgeom->SetGlobalID(index);
            m_hexGeoms[index] = hexgeom;
            return hexgeom;
        }

        // \brief Read the expansions given the XML file path.
        void MeshGraph::ReadExpansions(const std::string& infilename)
        {
          TiXmlDocument doc(infilename);
          bool loadOkay = doc.LoadFile();

          std::stringstream errstr;
          errstr << "Unable to load file: " << infilename << std::endl;
          errstr << "Reason: " << doc.ErrorDesc() << std::endl;
          errstr << "Position: Line " << doc.ErrorRow() << ", Column " << doc.ErrorCol() << std::endl;
          ASSERTL0(loadOkay, errstr.str());

          ReadExpansions(doc);
        }

        LibUtilities::BasisKeyVector DefineBasisKeyFromExpansionType(GeometrySharedPtr in,
                                     ExpansionType type,
                                     const int order)
        {
            LibUtilities::BasisKeyVector returnval;
            
            GeomShapeType shape= in->GetGeomShapeType();
            
            switch(type)
            {
					
            case eModified:
			{
                switch (shape)
                {
                case eSegment:
                    {
                        const LibUtilities::PointsKey pkey(order+1,LibUtilities::eGaussLobattoLegendre);
                        LibUtilities::BasisKey bkey(LibUtilities::eModified_A,order,pkey);
                        returnval.push_back(bkey);
                    }
                    break;
                case eQuadrilateral:
                    {
                        const LibUtilities::PointsKey pkey(order+1,LibUtilities::eGaussLobattoLegendre);
                        LibUtilities::BasisKey bkey(LibUtilities::eModified_A,order,pkey);
                        returnval.push_back(bkey);
                        returnval.push_back(bkey);
                    }
                    break;
                case eHexahedron:
                    {
                        const LibUtilities::PointsKey pkey(order+1,LibUtilities::eGaussLobattoLegendre);
                        LibUtilities::BasisKey bkey(LibUtilities::eModified_A,order,pkey);
                        returnval.push_back(bkey);
                        returnval.push_back(bkey);
                        returnval.push_back(bkey);
                    }
                    break;
                case eTriangle:
                    {
                        const LibUtilities::PointsKey pkey(order+1,LibUtilities::eGaussLobattoLegendre);
                        LibUtilities::BasisKey bkey(LibUtilities::eModified_A,order,pkey);
                        returnval.push_back(bkey);
                        
                        const LibUtilities::PointsKey pkey1(order,LibUtilities::eGaussRadauMAlpha1Beta0);
                        LibUtilities::BasisKey bkey1(LibUtilities::eModified_B,order,pkey1);
                        
                        returnval.push_back(bkey1);
                    }
                    break;
                case eTetrahedron:
                    {
                        const LibUtilities::PointsKey pkey(order+1,LibUtilities::eGaussLobattoLegendre);
                        LibUtilities::BasisKey bkey(LibUtilities::eModified_A,order,pkey);
                        returnval.push_back(bkey);
                        
                        const LibUtilities::PointsKey pkey1(order,LibUtilities::eGaussRadauMAlpha1Beta0);
                        LibUtilities::BasisKey bkey1(LibUtilities::eModified_B,order,pkey1);
                        returnval.push_back(bkey1);
                        
                        const LibUtilities::PointsKey pkey2(order,LibUtilities::eGaussRadauMAlpha2Beta0);
                        LibUtilities::BasisKey bkey2(LibUtilities::eModified_C,order,pkey2);
                        returnval.push_back(bkey2);
                    }
                    break;
                case ePrism:
                    {
                        const LibUtilities::PointsKey pkey(order+1,LibUtilities::eGaussLobattoLegendre);
                        LibUtilities::BasisKey bkey(LibUtilities::eModified_A,order,pkey);
                        returnval.push_back(bkey);
                        returnval.push_back(bkey);
                        
                        const LibUtilities::PointsKey pkey1(order,LibUtilities::eGaussRadauMAlpha1Beta0);
                        LibUtilities::BasisKey bkey1(LibUtilities::eModified_B,order,pkey1);
                        returnval.push_back(bkey1);
                        
                    }
                    break;
                default:
				{
                    ASSERTL0(false,"Expansion not defined in switch  for this shape");
				}
				break;
                }
			}
			break;
            
			case eGLL_Lagrange:
            {
                
                switch(shape)
                {
					case eSegment:
                    {
                        const LibUtilities::PointsKey pkey(order+1,LibUtilities::eGaussLobattoLegendre);
                        LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange,order,pkey);
                        returnval.push_back(bkey);
					}
                    break;
                    case eQuadrilateral:
					{
						const LibUtilities::PointsKey pkey(order+1,LibUtilities::eGaussLobattoLegendre);
						LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange,order,pkey);
						returnval.push_back(bkey);
						returnval.push_back(bkey);
					}
					break;
					case eTriangle: // define with corrects points key
						            // and change to Ortho on
                                    // construction
					{
						const LibUtilities::PointsKey pkey(order+1,LibUtilities::eGaussLobattoLegendre);
						LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange,order,pkey);
						returnval.push_back(bkey);
                    
						const LibUtilities::PointsKey pkey1(order,LibUtilities::eGaussRadauMAlpha1Beta0);
						LibUtilities::BasisKey bkey1(LibUtilities::eOrtho_B,order,pkey1);
						returnval.push_back(bkey1);
					}
					break;
                    default:
					{
						ASSERTL0(false,"Expansion not defined in switch  for this shape");
					}
                    break;
                }
            }
            break;
            
			case eOrthogonal:
			{
                switch (shape)
                {
                case eSegment:
				{
					const LibUtilities::PointsKey pkey(order+1,LibUtilities::eGaussLobattoLegendre);
					LibUtilities::BasisKey bkey(LibUtilities::eOrtho_A,order,pkey);
                        
					returnval.push_back(bkey);
				}
				break;
                case eTriangle:
                {
                    const LibUtilities::PointsKey pkey(order+1,LibUtilities::eGaussLobattoLegendre);
                    LibUtilities::BasisKey bkey(LibUtilities::eOrtho_A,order,pkey);
                    
                    returnval.push_back(bkey);
                    
                    const LibUtilities::PointsKey pkey1(order,LibUtilities::eGaussRadauMAlpha1Beta0);
                    LibUtilities::BasisKey bkey1(LibUtilities::eOrtho_B,order,pkey1);
                    
                    returnval.push_back(bkey1);
                }
                break;
                case eQuadrilateral:
                {
                    const LibUtilities::PointsKey pkey(order+1,LibUtilities::eGaussLobattoLegendre);
                    LibUtilities::BasisKey bkey(LibUtilities::eOrtho_A,order,pkey);
                    
                    returnval.push_back(bkey);
                    returnval.push_back(bkey);
                }
                break;
                case eTetrahedron:
                {
                    const LibUtilities::PointsKey pkey(order+1,LibUtilities::eGaussLobattoLegendre);
                    LibUtilities::BasisKey bkey(LibUtilities::eOrtho_A,order,pkey);
                    
                    returnval.push_back(bkey);
                    
                    const LibUtilities::PointsKey pkey1(order,LibUtilities::eGaussRadauMAlpha1Beta0);
                    LibUtilities::BasisKey bkey1(LibUtilities::eOrtho_B,order,pkey1);
                    
                    returnval.push_back(bkey1);
                    
                    const LibUtilities::PointsKey pkey2(order,LibUtilities::eGaussRadauMAlpha2Beta0);
                    LibUtilities::BasisKey bkey2(LibUtilities::eOrtho_C,order,pkey2);
                }
                break;
                default:
				{
                    ASSERTL0(false,"Expansion not defined in switch  for this shape");