<commit_before>////////////////////////////////////////////////////////////////////////////////
//
//  File:  $Source: /usr/sci/projects/Nektar/cvs/Nektar++/library/SpatialDomains/MeshGraph2D.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 <SpatialDomains/MeshGraph2D.h>
#include <SpatialDomains/SegGeom.h>
#include <SpatialDomains/ParseUtils.hpp>

namespace Nektar
{
    namespace SpatialDomains
    {
        MeshGraph2D::MeshGraph2D()
        {
        }

        MeshGraph2D::~MeshGraph2D()
        {
        }

        void MeshGraph2D::ReadGeometry(std::string &infilename)
        {
            TiXmlDocument doc(infilename);
            bool loadOkay = doc.LoadFile();

            ASSERTL0(loadOkay, (std::string("Unable to load file: ") + 
                                infilename).c_str());

            ReadGeometry(doc);
        }

        // \brief Read segments (and general MeshGraph) given TiXmlDocument.
        void MeshGraph2D::ReadGeometry(TiXmlDocument &doc)
        {
            // Read mesh first
            MeshGraph::ReadGeometry(doc);
            TiXmlHandle docHandle(&doc);

            TiXmlNode* node = NULL;
            TiXmlElement* mesh = NULL;

            /// Look for all geometry related data in GEOMETRY block.
            mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();

            ASSERTL0(mesh, "Unable to find GEOMETRY tag in file.");

            ReadCurves(doc);
            ReadEdges(doc);
            ReadElements(doc);
            ReadComposites(doc);
            ReadDomain(doc);
        }

        void MeshGraph2D::ReadEdges(TiXmlDocument &doc)
        {
            /// We know we have it since we made it this far.
            TiXmlHandle docHandle(&doc);
            TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
            TiXmlElement* field = NULL;

            /// Look for elements in ELEMENT block.
            field = mesh->FirstChildElement("EDGE");

            ASSERTL0(field, "Unable to find EDGE tag in file.");

            /// All elements are of the form: "<E ID="#"> ... </E>", with
            /// ? being the element type.
            /// Read the ID field first.
            TiXmlElement *edge = field->FirstChildElement("E");

            /// Since all edge data is one big text block, we need to
            /// accumulate all TEXT data and then parse it.  This
            /// approach effectively skips all comments or other node
            /// types since we only care about the edge list.  We
            /// cannot handle missing edge numbers as we could with
            /// missing element numbers due to the text block format.
            std::string edgeStr;
            int i,indx;
            int err = 0;
            int nextEdgeNumber = -1;

            // Curved Edges
            map<int, int> edge_curved;
            for(i = 0; i < m_curvededges.size(); ++i)
            {
                edge_curved[m_curvededges[i]->m_curveID] = i;
            }

            while(edge)
            {
                nextEdgeNumber++;

                int err = edge->QueryIntAttribute("ID",&indx);
                ASSERTL0(err == TIXML_SUCCESS, "Unable to read edge attribute ID.");
                ASSERTL0(indx == nextEdgeNumber, "Edge IDs must begin with zero and be sequential.");

                TiXmlNode *child = edge->FirstChild();
                edgeStr.clear();
                if (child->Type() == TiXmlNode::TEXT)
                {
                    edgeStr += child->ToText()->ValueStr();
                }

                /// Now parse out the edges, three fields at a time.
                int vertex1, vertex2;
                std::istringstream edgeDataStrm(edgeStr.c_str());

                try
                {
                    while (!edgeDataStrm.fail())
                    {
                        edgeDataStrm >> vertex1 >> vertex2;

                        // Must check after the read because we may be
                        // at the end and not know it.  If we are at
                        // the end we will add a duplicate of the last
                        // entry if we don't check here.
                        if (!edgeDataStrm.fail())
                        {
                            VertexComponentSharedPtr vertices[2] = {GetVertex(vertex1), GetVertex(vertex2)};
                            
                            SegGeomSharedPtr edge; 

                            if(edge_curved.count(indx) == 0)
                            {
                                edge = MemoryManager<SegGeom>::AllocateSharedPtr(indx, m_SpaceDimension, vertices);
                            }
                            else
                            {
                                edge = MemoryManager<SegGeom>::AllocateSharedPtr(indx, m_SpaceDimension, vertices, m_curvededges[edge_curved.find(indx)->second]);
                            }
                            
                            m_seggeoms.push_back(edge);
                        }
                    }
                }
                catch(...)
                {
                    NEKERROR(ErrorUtil::efatal, (std::string("Unable to read edge data: ") + edgeStr).c_str());
                }

                edge = edge->NextSiblingElement("E");
            }

        }

        void MeshGraph2D::ReadElements(TiXmlDocument &doc)
        {
            /// We know we have it since we made it this far.
            TiXmlHandle docHandle(&doc);
            TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
            TiXmlElement* field = NULL;

            /// Look for elements in ELEMENT block.
            field = mesh->FirstChildElement("ELEMENT");

            ASSERTL0(field, "Unable to find ELEMENT tag in file.");

            int nextElementNumber = -1;

            /// All elements are of the form: "<? ID="#"> ... </?>", with
            /// ? being the element type.

            TiXmlElement *element = field->FirstChildElement();

            while (element)
            {
                std::string elementType(element->ValueStr());

                ASSERTL0(elementType == "Q" || elementType == "T",
                    (std::string("Unknown 2D element type: ") + elementType).c_str());

                /// These should be ordered.
                nextElementNumber++;

                /// Read id attribute.
                int indx;
                int err = element->QueryIntAttribute("ID", &indx);
                ASSERTL0(err == TIXML_SUCCESS, "Unable to read element attribute ID.");
                ASSERTL0(indx == nextElementNumber, "Element IDs must begin with zero and be sequential.");

                /// Read text element description.
                TiXmlNode* elementChild = element->FirstChild();
                std::string elementStr;
                while(elementChild)
                {
                    if (elementChild->Type() == TiXmlNode::TEXT)
                    {
                        elementStr += elementChild->ToText()->ValueStr();
                    }
                    elementChild = elementChild->NextSibling();
                }

                ASSERTL0(!elementStr.empty(), "Unable to read element description body.");

                /// Parse out the element components corresponding to type of element.
                if (elementType == "T")
                {
                    // Read three edge numbers
                    int edge1, edge2, edge3;
                    std::istringstream elementDataStrm(elementStr.c_str());

                    try
                    {
                        elementDataStrm >> edge1;
                        elementDataStrm >> edge2;
                        elementDataStrm >> edge3;

                        ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for TRIANGLE: ") + elementStr).c_str());

                        /// Create a TriGeom to hold the new definition.
                        SegGeomSharedPtr edges[TriGeom::kNedges] = 
                        {
                            GetSegGeom(edge1),
                            GetSegGeom(edge2),
                            GetSegGeom(edge3)
                        };

                        StdRegions::EdgeOrientation edgeorient[TriGeom::kNedges] = 
                        {
                            SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
                            SegGeom::GetEdgeOrientation(*edges[1], *edges[2]), 
                            SegGeom::GetEdgeOrientation(*edges[2], *edges[0])
                        };

                        TriGeomSharedPtr trigeom(MemoryManager<TriGeom>::AllocateSharedPtr(indx,edges,edgeorient));
                        trigeom->SetGlobalID(indx);

                        m_trigeoms.push_back(trigeom);
                    }
                    catch(...)
                    {
                        NEKERROR(ErrorUtil::efatal,
                            (std::string("Unable to read element data for TRIANGLE: ") + elementStr).c_str());
                    }
                }
                else if (elementType == "Q")
                {
                    // Read four edge numbers
                    int edge1, edge2, edge3, edge4;
                    std::istringstream elementDataStrm(elementStr.c_str());

                    try
                    {
                        elementDataStrm >> edge1;
                        elementDataStrm >> edge2;
                        elementDataStrm >> edge3;
                        elementDataStrm >> edge4;

                        ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for QUAD: ") + elementStr).c_str());

                        /// Create a QuadGeom to hold the new definition.
                        SegGeomSharedPtr edges[QuadGeom::kNedges] = 
                        {GetSegGeom(edge1),GetSegGeom(edge2),
                         GetSegGeom(edge3),GetSegGeom(edge4)};

                        StdRegions::EdgeOrientation edgeorient[QuadGeom::kNedges] =
                        {
                            SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
                            SegGeom::GetEdgeOrientation(*edges[1], *edges[2]),
                            SegGeom::GetEdgeOrientation(*edges[2], *edges[3]),
                            SegGeom::GetEdgeOrientation(*edges[3], *edges[0])
                        };

                        //QuadGeomSharedPtr quadgeom(new QuadGeom(edges, edgeorient));
                        QuadGeomSharedPtr quadgeom  = MemoryManager<QuadGeom>::AllocateSharedPtr(indx, edges, edgeorient);
                        quadgeom->SetGlobalID(indx);
                        
                        m_quadgeoms.push_back(quadgeom);

                    }
                    catch(...)
                    {
                        NEKERROR(ErrorUtil::efatal,(std::string("Unable to read element data for QUAD: ") + elementStr).c_str());
                    }
                }

                /// Keep looking
                element = element->NextSiblingElement();
            }
        }

        void MeshGraph2D::ReadComposites(TiXmlDocument &doc)
        {
            TiXmlHandle docHandle(&doc);

            /// We know we have it since we made it this far.
            TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
            TiXmlElement* field = NULL;

            ASSERTL0(mesh, "Unable to find GEOMETRY tag in file.");

            /// Look for elements in ELEMENT block.
            field = mesh->FirstChildElement("COMPOSITE");

            ASSERTL0(field, "Unable to find COMPOSITE tag in file.");

            int nextCompositeNumber = -1;

            /// All elements are of the form: "<C ID = "N"> ... </C>".

            /// Read the ID field first.
            TiXmlElement *composite = field->FirstChildElement("C");

            while (composite)
            {
                nextCompositeNumber++;

                int indx;
                int err = composite->QueryIntAttribute("ID", &indx);
                ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
                ASSERTL0(indx == nextCompositeNumber, "Composite IDs must begin with zero and be sequential.");

                TiXmlNode* compositeChild = composite->FirstChild();
                // This is primarily to skip comments that may be present.
                // Comments appear as nodes just like elements.
                // We are specifically looking for text in the body
                // of the definition.
                while(compositeChild && compositeChild->Type() != TiXmlNode::TEXT)
                {
                    compositeChild = compositeChild->NextSibling();
                }

                ASSERTL0(compositeChild, "Unable to read composite definition body.");
                std::string compositeStr = compositeChild->ToText()->ValueStr();

                /// Parse out the element components corresponding to type of element.

                std::istringstream compositeDataStrm(compositeStr.c_str());

                try
                {
                    bool first = true;
                    std::string prevCompositeElementStr;

                    while (!compositeDataStrm.fail())
                    {
                        std::string compositeElementStr;
                        compositeDataStrm >> compositeElementStr;

                        if (!compositeDataStrm.fail())
                        {
                            if (first)
                            {
                                first = false;

                                Composite curVector(MemoryManager<GeometryVector>::AllocateSharedPtr());
                                m_MeshCompositeVector.push_back(curVector);
                            }

                            if (compositeElementStr.length() > 0)
                            {
                                ResolveGeomRef(prevCompositeElementStr, compositeElementStr);
                            }
                            prevCompositeElementStr = compositeElementStr;
                        }
                    }
                }
                catch(...)
                {
                    NEKERROR(ErrorUtil::efatal,
                        (std::string("Unable to read COMPOSITE data for composite: ") + compositeStr).c_str());
                }

                /// Keep looking
                composite = composite->NextSiblingElement("C");
            }
        }


        SegGeomSharedPtr MeshGraph2D::GetSegGeom(int eID)
        {
            SegGeomSharedPtr returnval;

            if (eID >= 0 && eID < int(m_seggeoms.size()))
            {
                returnval = m_seggeoms[eID];
            }

            return returnval;
        };


        // Take the string that is the composite reference and find the
        // pointer to the Geometry object corresponding to it.

        // The only allowable combinations of previous and current items
        // are V (0D); E (1D); and T and Q (2D).  Only elements of the same
        // dimension are allowed to be grouped.
        void MeshGraph2D::ResolveGeomRef(const std::string &prevToken, const std::string &token)
        {
            try
            {
                std::istringstream tokenStream(token);
                std::istringstream prevTokenStream(prevToken);

                char type;
                char prevType;

                tokenStream >> type;

                std::string::size_type indxBeg = token.find_first_of('[') + 1;
                std::string::size_type indxEnd = token.find_last_of(']') - 1;

                ASSERTL0(indxBeg <= indxEnd, (std::string("Error reading index definition:") + token).c_str());

                std::string indxStr = token.substr(indxBeg, indxEnd - indxBeg + 1);

                std::vector<unsigned int> seqVector;
                std::vector<unsigned int>::iterator seqIter;

                bool err = ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector);

                ASSERTL0(err, (std::string("Error reading composite elements: ") + indxStr).c_str());

                prevTokenStream >> prevType;

                // All composites must be of the same dimension.
                bool validSequence = (prevToken.empty() ||         // No previous, then current is just fine.
                    (type == 'V' && prevType == 'V') ||
                    (type == 'E' && prevType == 'E') ||
                    ((type == 'T' || type == 'Q') &&
                    (prevType == 'T' || prevType == 'Q')));

                ASSERTL0(validSequence, (std::string("Invalid combination of composite items: ")
                    + type + " and " + prevType + ".").c_str()); 

                switch(type)
                {
                case 'E':   // Edge
                    for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
                    {
                        if (*seqIter >= m_seggeoms.size())
                        {
                            char errStr[16] = "";
                            ::sprintf(errStr, "%d", *seqIter);
                            NEKERROR(ErrorUtil::ewarning, (std::string("Unknown edge index: ") + errStr).c_str());
                        }
                        else
                        {
                            m_MeshCompositeVector.back()->push_back(m_seggeoms[*seqIter]);
                        }
                    }
                    break;

                case 'T':   // Triangle
                    for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
                    {
                        if (*seqIter >= m_trigeoms.size())
                        {
                            char errStr[16] = "";
                            ::sprintf(errStr, "%d", *seqIter);
                            NEKERROR(ErrorUtil::ewarning, (std::string("Unknown triangle index: ") + errStr).c_str());
                        }
                        else
                        {
                            m_MeshCompositeVector.back()->push_back(m_trigeoms[*seqIter]);
                        }
                    }
                    break;

                case 'Q':   // Quad
                    for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
                    {
                        if (*seqIter >= m_quadgeoms.size())
                        {
                            char errStr[16] = "";
                            ::sprintf(errStr, "%d", *seqIter);
                            NEKERROR(ErrorUtil::ewarning, (std::string("Unknown quad index: ") + errStr).c_str());
                        }
                        else
                        {
                            m_MeshCompositeVector.back()->push_back(m_quadgeoms[*seqIter]);
                        }
                    }
                    break;

                case 'V':   // Vertex
                    for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
                    {
                        if (*seqIter >= m_vertset.size())
                        {
                            char errStr[16] = "";
                            ::sprintf(errStr, "%d", *seqIter);
                            NEKERROR(ErrorUtil::ewarning, (std::string("Unknown vertex index: ") + errStr).c_str());
                        }
                        else
                        {
                            m_MeshCompositeVector.back()->push_back(m_vertset[*seqIter]);
                        }
                    }
                    break;

                default:
                    NEKERROR(ErrorUtil::efatal, (std::string("Unrecognized composite token: ") + token).c_str());
                }
            }
            catch(...)
            {
                NEKERROR(ErrorUtil::efatal, (std::string("Problem processing composite token: ") + token).c_str());
            }

            return;
        }
        
        ElementEdgeVectorSharedPtr MeshGraph2D::GetElementsFromEdge(Geometry1DSharedPtr edge)
        {
            SegGeomSharedPtr Sedge;

            if(!(Sedge = boost::dynamic_pointer_cast<SegGeom>(edge)))
            {
                ASSERTL0(false,"Dynamics cast failed");

            }
            return GetElementsFromEdge(Sedge);
            
        }
        ElementEdgeVectorSharedPtr MeshGraph2D::GetElementsFromEdge(SegGeomSharedPtr edge)
        {
            // Search tris and quads
            // Need to iterate through vectors because there may be multiple
            // occurrences.
            ElementEdgeSharedPtr elementEdge;
            //TriGeomVectorIter triIter;

            ElementEdgeVectorSharedPtr returnval = MemoryManager<ElementEdgeVector>::AllocateSharedPtr();

            CompositeVectorIter compIter;
            TriGeomSharedPtr triGeomShPtr;
            QuadGeomSharedPtr quadGeomShPtr;

            GeometryVectorIter geomIter;

            for (compIter = m_Domain.begin(); compIter != m_Domain.end(); ++compIter)
            {
                for (geomIter = (*compIter)->begin(); geomIter != (*compIter)->end(); ++geomIter)
                {
                    triGeomShPtr = boost::dynamic_pointer_cast<TriGeom>(*geomIter);
                    quadGeomShPtr = boost::dynamic_pointer_cast<QuadGeom>(*geomIter);

                    if (triGeomShPtr || quadGeomShPtr)
                    {
                        int edgeNum;                        
                        if (triGeomShPtr)
                        {
                            if ((edgeNum = triGeomShPtr->WhichEdge(edge)) > -1)
                            {
                                elementEdge = MemoryManager<ElementEdge>::AllocateSharedPtr();
                                elementEdge->m_Element = triGeomShPtr;
                                elementEdge->m_EdgeIndx = edgeNum;
                                returnval->push_back(elementEdge);
                            }
                        }
                        else if (quadGeomShPtr)
                        {
                            if ((edgeNum = quadGeomShPtr->WhichEdge(edge)) > -1)
                            {
                                elementEdge = MemoryManager<ElementEdge>::AllocateSharedPtr();
                                elementEdge->m_Element = quadGeomShPtr;
                                elementEdge->m_EdgeIndx = edgeNum;
                                returnval->push_back(elementEdge);
                            }
                        }
                    }
                }         
            }

            //for(triIter = m_trigeoms.begin(); triIter != m_trigeoms.end(); ++triIter)
            //{
            //    int edgeNum;
            //    if ((edgeNum = (*triIter)->WhichEdge(edge)) > -1)
            //    {
            //        elementEdge = MemoryManager<ElementEdge>::AllocateSharedPtr();
            //        elementEdge->m_Element = *triIter;
            //        elementEdge->m_EdgeIndx = edgeNum;
            //        returnval->push_back(elementEdge);
            //    }
            //}

            //QuadGeomVector::iterator quadIter;

            //for(quadIter = m_quadgeoms.begin(); quadIter != m_quadgeoms.end(); ++quadIter)
            //{
            //    int edgeNum;
            //    if ((edgeNum = (*quadIter)->WhichEdge(edge)) > -1)
            //    {
            //        elementEdge = MemoryManager<ElementEdge>::AllocateSharedPtr();
            //        elementEdge->m_Element = *quadIter;
            //        elementEdge->m_EdgeIndx = edgeNum;
            //        returnval->push_back(elementEdge);
            //    }
            //}

            return returnval;
        }

        LibUtilities::BasisKey MeshGraph2D::GetEdgeBasisKey(SegGeomSharedPtr edge)
        {
            ElementEdgeVectorSharedPtr elements = GetElementsFromEdge(edge);
            // Perhaps, a check should be done here to ensure that in case 
            // elements->size!=1, all elements to which the edge belongs have the same type
            // and order of expansion such that no confusion can arise.
            ExpansionShPtr expansion = GetExpansion((*elements)[0]->m_Element);
            

            int edge_id = (*elements)[0]->m_EdgeIndx;

            if((*elements)[0]->m_Element->GetGeomShapeType() == eTriangle)
            {
                edge_id = (edge_id)? 1:0;
            }
            else
            {
                edge_id = edge_id%2;
            }

            int nummodes  = expansion->m_BasisKeyVector[edge_id].GetNumModes();
            int numpoints = expansion->m_BasisKeyVector[edge_id].GetNumPoints();
            
	    if((*elements)[0]->m_Element->GetGeomShapeType() == eTriangle)
	      {
		// Use edge 0 to define basis of order relevant to edge 
                switch(expansion->m_BasisKeyVector[edge_id].GetBasisType())
		  {
		  case LibUtilities::eModified_B:
                    {
		      switch(expansion->m_BasisKeyVector[edge_id].GetPointsType())
                        {
                        case LibUtilities::eGaussRadauMAlpha1Beta0:
			  {                            
			    const LibUtilities::PointsKey pkey(numpoints+1,LibUtilities::eGaussLobattoLegendre);
			    return LibUtilities::BasisKey(expansion->m_BasisKeyVector[0].GetBasisType(),nummodes,pkey);
			  }
			  break;
			default:
			  ASSERTL0(false,"Unexpected points distribution");
			  break;
			}
		    }
		    break;
		  case LibUtilities::eModified_A:
                    {
		      switch(expansion->m_BasisKeyVector[edge_id].GetPointsType())
                        {
                        case LibUtilities::eGaussLobattoLegendre:
			  {                            
			    const LibUtilities::PointsKey pkey(numpoints,LibUtilities::eGaussLobattoLegendre);
			    return LibUtilities::BasisKey(expansion->m_BasisKeyVector[0].GetBasisType(),nummodes,pkey);
			  }
			  break;
			default:
			  ASSERTL0(false,"Unexpected points distribution");
			  break;
			}
		    }
		    break;
		  default:
		    ASSERTL0(false,"Unexpected basis distribution");
		  }
	      }
	    else
	      {
		// Quadrilateral
		const LibUtilities::PointsKey pkey(numpoints,expansion->m_BasisKeyVector[edge_id].GetPointsType());
		return LibUtilities::BasisKey(expansion->m_BasisKeyVector[edge_id].GetBasisType(),nummodes,pkey);
	      }
        }
      
    }; //end of namespace
}; //end of namespace

//
// $Log: MeshGraph2D.cpp,v $
// Revision 1.39  2009/07/02 13:32:24  sehunchun
// *** empty log message ***
//
// Revision 1.38  2009/07/02 11:39:44  sehunchun
// Modification for 2D gemoetry embedded in 3D
//<commit_msg>Fixed visual studio compiler warning.<commit_after>////////////////////////////////////////////////////////////////////////////////
//
//  File:  $Source: /usr/sci/projects/Nektar/cvs/Nektar++/library/SpatialDomains/MeshGraph2D.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 <SpatialDomains/MeshGraph2D.h>
#include <SpatialDomains/SegGeom.h>
#include <SpatialDomains/ParseUtils.hpp>

namespace Nektar
{
    namespace SpatialDomains
    {
        MeshGraph2D::MeshGraph2D()
        {
        }

        MeshGraph2D::~MeshGraph2D()
        {
        }

        void MeshGraph2D::ReadGeometry(std::string &infilename)
        {
            TiXmlDocument doc(infilename);
            bool loadOkay = doc.LoadFile();

            ASSERTL0(loadOkay, (std::string("Unable to load file: ") + 
                                infilename).c_str());

            ReadGeometry(doc);
        }

        // \brief Read segments (and general MeshGraph) given TiXmlDocument.
        void MeshGraph2D::ReadGeometry(TiXmlDocument &doc)
        {
            // Read mesh first
            MeshGraph::ReadGeometry(doc);
            TiXmlHandle docHandle(&doc);

            TiXmlNode* node = NULL;
            TiXmlElement* mesh = NULL;

            /// Look for all geometry related data in GEOMETRY block.
            mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();

            ASSERTL0(mesh, "Unable to find GEOMETRY tag in file.");

            ReadCurves(doc);
            ReadEdges(doc);
            ReadElements(doc);
            ReadComposites(doc);
            ReadDomain(doc);
        }

        void MeshGraph2D::ReadEdges(TiXmlDocument &doc)
        {
            /// We know we have it since we made it this far.
            TiXmlHandle docHandle(&doc);
            TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
            TiXmlElement* field = NULL;

            /// Look for elements in ELEMENT block.
            field = mesh->FirstChildElement("EDGE");

            ASSERTL0(field, "Unable to find EDGE tag in file.");

            /// All elements are of the form: "<E ID="#"> ... </E>", with
            /// ? being the element type.
            /// Read the ID field first.
            TiXmlElement *edge = field->FirstChildElement("E");

            /// Since all edge data is one big text block, we need to
            /// accumulate all TEXT data and then parse it.  This
            /// approach effectively skips all comments or other node
            /// types since we only care about the edge list.  We
            /// cannot handle missing edge numbers as we could with
            /// missing element numbers due to the text block format.
            std::string edgeStr;
            int i,indx;
            int err = 0;
            int nextEdgeNumber = -1;

            // Curved Edges
            map<int, int> edge_curved;
            for(i = 0; i < m_curvededges.size(); ++i)
            {
                edge_curved[m_curvededges[i]->m_curveID] = i;
            }

            while(edge)
            {
                nextEdgeNumber++;

                int err = edge->QueryIntAttribute("ID",&indx);
                ASSERTL0(err == TIXML_SUCCESS, "Unable to read edge attribute ID.");
                ASSERTL0(indx == nextEdgeNumber, "Edge IDs must begin with zero and be sequential.");

                TiXmlNode *child = edge->FirstChild();
                edgeStr.clear();
                if (child->Type() == TiXmlNode::TEXT)
                {
                    edgeStr += child->ToText()->ValueStr();
                }

                /// Now parse out the edges, three fields at a time.
                int vertex1, vertex2;
                std::istringstream edgeDataStrm(edgeStr.c_str());

                try
                {
                    while (!edgeDataStrm.fail())
                    {
                        edgeDataStrm >> vertex1 >> vertex2;

                        // Must check after the read because we may be
                        // at the end and not know it.  If we are at
                        // the end we will add a duplicate of the last
                        // entry if we don't check here.
                        if (!edgeDataStrm.fail())
                        {
                            VertexComponentSharedPtr vertices[2] = {GetVertex(vertex1), GetVertex(vertex2)};
                            
                            SegGeomSharedPtr edge; 

                            if(edge_curved.count(indx) == 0)
                            {
                                edge = MemoryManager<SegGeom>::AllocateSharedPtr(indx, m_SpaceDimension, vertices);
                            }
                            else
                            {
                                edge = MemoryManager<SegGeom>::AllocateSharedPtr(indx, m_SpaceDimension, vertices, m_curvededges[edge_curved.find(indx)->second]);
                            }
                            
                            m_seggeoms.push_back(edge);
                        }
                    }
                }
                catch(...)
                {
                    NEKERROR(ErrorUtil::efatal, (std::string("Unable to read edge data: ") + edgeStr).c_str());
                }

                edge = edge->NextSiblingElement("E");
            }

        }

        void MeshGraph2D::ReadElements(TiXmlDocument &doc)
        {
            /// We know we have it since we made it this far.
            TiXmlHandle docHandle(&doc);
            TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
            TiXmlElement* field = NULL;

            /// Look for elements in ELEMENT block.
            field = mesh->FirstChildElement("ELEMENT");

            ASSERTL0(field, "Unable to find ELEMENT tag in file.");

            int nextElementNumber = -1;

            /// All elements are of the form: "<? ID="#"> ... </?>", with
            /// ? being the element type.

            TiXmlElement *element = field->FirstChildElement();

            while (element)
            {
                std::string elementType(element->ValueStr());

                ASSERTL0(elementType == "Q" || elementType == "T",
                    (std::string("Unknown 2D element type: ") + elementType).c_str());

                /// These should be ordered.
                nextElementNumber++;

                /// Read id attribute.
                int indx;
                int err = element->QueryIntAttribute("ID", &indx);
                ASSERTL0(err == TIXML_SUCCESS, "Unable to read element attribute ID.");
                ASSERTL0(indx == nextElementNumber, "Element IDs must begin with zero and be sequential.");

                /// Read text element description.
                TiXmlNode* elementChild = element->FirstChild();
                std::string elementStr;
                while(elementChild)
                {
                    if (elementChild->Type() == TiXmlNode::TEXT)
                    {
                        elementStr += elementChild->ToText()->ValueStr();
                    }
                    elementChild = elementChild->NextSibling();
                }

                ASSERTL0(!elementStr.empty(), "Unable to read element description body.");

                /// Parse out the element components corresponding to type of element.
                if (elementType == "T")
                {
                    // Read three edge numbers
                    int edge1, edge2, edge3;
                    std::istringstream elementDataStrm(elementStr.c_str());

                    try
                    {
                        elementDataStrm >> edge1;
                        elementDataStrm >> edge2;
                        elementDataStrm >> edge3;

                        ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for TRIANGLE: ") + elementStr).c_str());

                        /// Create a TriGeom to hold the new definition.
                        SegGeomSharedPtr edges[TriGeom::kNedges] = 
                        {
                            GetSegGeom(edge1),
                            GetSegGeom(edge2),
                            GetSegGeom(edge3)
                        };

                        StdRegions::EdgeOrientation edgeorient[TriGeom::kNedges] = 
                        {
                            SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
                            SegGeom::GetEdgeOrientation(*edges[1], *edges[2]), 
                            SegGeom::GetEdgeOrientation(*edges[2], *edges[0])
                        };

                        TriGeomSharedPtr trigeom(MemoryManager<TriGeom>::AllocateSharedPtr(indx,edges,edgeorient));
                        trigeom->SetGlobalID(indx);

                        m_trigeoms.push_back(trigeom);
                    }
                    catch(...)
                    {
                        NEKERROR(ErrorUtil::efatal,
                            (std::string("Unable to read element data for TRIANGLE: ") + elementStr).c_str());
                    }
                }
                else if (elementType == "Q")
                {
                    // Read four edge numbers
                    int edge1, edge2, edge3, edge4;
                    std::istringstream elementDataStrm(elementStr.c_str());

                    try
                    {
                        elementDataStrm >> edge1;
                        elementDataStrm >> edge2;
                        elementDataStrm >> edge3;
                        elementDataStrm >> edge4;

                        ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for QUAD: ") + elementStr).c_str());

                        /// Create a QuadGeom to hold the new definition.
                        SegGeomSharedPtr edges[QuadGeom::kNedges] = 
                        {GetSegGeom(edge1),GetSegGeom(edge2),
                         GetSegGeom(edge3),GetSegGeom(edge4)};

                        StdRegions::EdgeOrientation edgeorient[QuadGeom::kNedges] =
                        {
                            SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
                            SegGeom::GetEdgeOrientation(*edges[1], *edges[2]),
                            SegGeom::GetEdgeOrientation(*edges[2], *edges[3]),
                            SegGeom::GetEdgeOrientation(*edges[3], *edges[0])
                        };

                        //QuadGeomSharedPtr quadgeom(new QuadGeom(edges, edgeorient));
                        QuadGeomSharedPtr quadgeom  = MemoryManager<QuadGeom>::AllocateSharedPtr(indx, edges, edgeorient);
                        quadgeom->SetGlobalID(indx);
                        
                        m_quadgeoms.push_back(quadgeom);

                    }
                    catch(...)
                    {
                        NEKERROR(ErrorUtil::efatal,(std::string("Unable to read element data for QUAD: ") + elementStr).c_str());
                    }
                }

                /// Keep looking
                element = element->NextSiblingElement();
            }
        }

        void MeshGraph2D::ReadComposites(TiXmlDocument &doc)
        {
            TiXmlHandle docHandle(&doc);

            /// We know we have it since we made it this far.
            TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
            TiXmlElement* field = NULL;

            ASSERTL0(mesh, "Unable to find GEOMETRY tag in file.");

            /// Look for elements in ELEMENT block.
            field = mesh->FirstChildElement("COMPOSITE");

            ASSERTL0(field, "Unable to find COMPOSITE tag in file.");

            int nextCompositeNumber = -1;

            /// All elements are of the form: "<C ID = "N"> ... </C>".

            /// Read the ID field first.
            TiXmlElement *composite = field->FirstChildElement("C");

            while (composite)
            {
                nextCompositeNumber++;

                int indx;
                int err = composite->QueryIntAttribute("ID", &indx);
                ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
                ASSERTL0(indx == nextCompositeNumber, "Composite IDs must begin with zero and be sequential.");

                TiXmlNode* compositeChild = composite->FirstChild();
                // This is primarily to skip comments that may be present.
                // Comments appear as nodes just like elements.
                // We are specifically looking for text in the body
                // of the definition.
                while(compositeChild && compositeChild->Type() != TiXmlNode::TEXT)
                {
                    compositeChild = compositeChild->NextSibling();
                }

                ASSERTL0(compositeChild, "Unable to read composite definition body.");
                std::string compositeStr = compositeChild->ToText()->ValueStr();

                /// Parse out the element components corresponding to type of element.

                std::istringstream compositeDataStrm(compositeStr.c_str());

                try
                {
                    bool first = true;
                    std::string prevCompositeElementStr;

                    while (!compositeDataStrm.fail())
                    {
                        std::string compositeElementStr;
                        compositeDataStrm >> compositeElementStr;

                        if (!compositeDataStrm.fail())
                        {
                            if (first)
                            {
                                first = false;

                                Composite curVector(MemoryManager<GeometryVector>::AllocateSharedPtr());
                                m_MeshCompositeVector.push_back(curVector);
                            }

                            if (compositeElementStr.length() > 0)
                            {
                                ResolveGeomRef(prevCompositeElementStr, compositeElementStr);
                            }
                            prevCompositeElementStr = compositeElementStr;
                        }
                    }
                }
                catch(...)
                {
                    NEKERROR(ErrorUtil::efatal,
                        (std::string("Unable to read COMPOSITE data for composite: ") + compositeStr).c_str());
                }

                /// Keep looking
                composite = composite->NextSiblingElement("C");
            }
        }


        SegGeomSharedPtr MeshGraph2D::GetSegGeom(int eID)
        {
            SegGeomSharedPtr returnval;

            if (eID >= 0 && eID < int(m_seggeoms.size()))
            {
                returnval = m_seggeoms[eID];
            }

            return returnval;
        };


        // Take the string that is the composite reference and find the
        // pointer to the Geometry object corresponding to it.

        // The only allowable combinations of previous and current items
        // are V (0D); E (1D); and T and Q (2D).  Only elements of the same
        // dimension are allowed to be grouped.
        void MeshGraph2D::ResolveGeomRef(const std::string &prevToken, const std::string &token)
        {
            try
            {
                std::istringstream tokenStream(token);
                std::istringstream prevTokenStream(prevToken);

                char type;
                char prevType;

                tokenStream >> type;

                std::string::size_type indxBeg = token.find_first_of('[') + 1;
                std::string::size_type indxEnd = token.find_last_of(']') - 1;

                ASSERTL0(indxBeg <= indxEnd, (std::string("Error reading index definition:") + token).c_str());

                std::string indxStr = token.substr(indxBeg, indxEnd - indxBeg + 1);

                std::vector<unsigned int> seqVector;
                std::vector<unsigned int>::iterator seqIter;

                bool err = ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector);

                ASSERTL0(err, (std::string("Error reading composite elements: ") + indxStr).c_str());

                prevTokenStream >> prevType;

                // All composites must be of the same dimension.
                bool validSequence = (prevToken.empty() ||         // No previous, then current is just fine.
                    (type == 'V' && prevType == 'V') ||
                    (type == 'E' && prevType == 'E') ||
                    ((type == 'T' || type == 'Q') &&
                    (prevType == 'T' || prevType == 'Q')));

                ASSERTL0(validSequence, (std::string("Invalid combination of composite items: ")
                    + type + " and " + prevType + ".").c_str()); 

                switch(type)
                {
                case 'E':   // Edge
                    for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
                    {
                        if (*seqIter >= m_seggeoms.size())
                        {
                            char errStr[16] = "";
                            ::sprintf(errStr, "%d", *seqIter);
                            NEKERROR(ErrorUtil::ewarning, (std::string("Unknown edge index: ") + errStr).c_str());
                        }
                        else
                        {
                            m_MeshCompositeVector.back()->push_back(m_seggeoms[*seqIter]);
                        }
                    }
                    break;

                case 'T':   // Triangle
                    for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
                    {
                        if (*seqIter >= m_trigeoms.size())
                        {
                            char errStr[16] = "";
                            ::sprintf(errStr, "%d", *seqIter);
                            NEKERROR(ErrorUtil::ewarning, (std::string("Unknown triangle index: ") + errStr).c_str());
                        }
                        else
                        {
                            m_MeshCompositeVector.back()->push_back(m_trigeoms[*seqIter]);
                        }
                    }
                    break;

                case 'Q':   // Quad
                    for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
                    {
                        if (*seqIter >= m_quadgeoms.size())
                        {
                            char errStr[16] = "";
                            ::sprintf(errStr, "%d", *seqIter);
                            NEKERROR(ErrorUtil::ewarning, (std::string("Unknown quad index: ") + errStr).c_str());
                        }
                        else
                        {
                            m_MeshCompositeVector.back()->push_back(m_quadgeoms[*seqIter]);
                        }
                    }
                    break;

                case 'V':   // Vertex
                    for (seqIter = seqVector.begin(); seqIter != seqVector.end(); ++seqIter)
                    {
                        if (*seqIter >= m_vertset.size())
                        {
                            char errStr[16] = "";
                            ::sprintf(errStr, "%d", *seqIter);
                            NEKERROR(ErrorUtil::ewarning, (std::string("Unknown vertex index: ") + errStr).c_str());
                        }
                        else
                        {
                            m_MeshCompositeVector.back()->push_back(m_vertset[*seqIter]);
                        }
                    }
                    break;

                default:
                    NEKERROR(ErrorUtil::efatal, (std::string("Unrecognized composite token: ") + token).c_str());
                }
            }
            catch(...)
            {
                NEKERROR(ErrorUtil::efatal, (std::string("Problem processing composite token: ") + token).c_str());
            }

            return;
        }
        
        ElementEdgeVectorSharedPtr MeshGraph2D::GetElementsFromEdge(Geometry1DSharedPtr edge)
        {
            SegGeomSharedPtr Sedge;

            if(!(Sedge = boost::dynamic_pointer_cast<SegGeom>(edge)))
            {
                ASSERTL0(false,"Dynamics cast failed");

            }
            return GetElementsFromEdge(Sedge);
            
        }
        ElementEdgeVectorSharedPtr MeshGraph2D::GetElementsFromEdge(SegGeomSharedPtr edge)
        {
            // Search tris and quads
            // Need to iterate through vectors because there may be multiple
            // occurrences.
            ElementEdgeSharedPtr elementEdge;
            //TriGeomVectorIter triIter;

            ElementEdgeVectorSharedPtr returnval = MemoryManager<ElementEdgeVector>::AllocateSharedPtr();

            CompositeVectorIter compIter;
            TriGeomSharedPtr triGeomShPtr;
            QuadGeomSharedPtr quadGeomShPtr;

            GeometryVectorIter geomIter;

            for (compIter = m_Domain.begin(); compIter != m_Domain.end(); ++compIter)
            {
                for (geomIter = (*compIter)->begin(); geomIter != (*compIter)->end(); ++geomIter)
                {
                    triGeomShPtr = boost::dynamic_pointer_cast<TriGeom>(*geomIter);
                    quadGeomShPtr = boost::dynamic_pointer_cast<QuadGeom>(*geomIter);

                    if (triGeomShPtr || quadGeomShPtr)
                    {
                        int edgeNum;                        
                        if (triGeomShPtr)
                        {
                            if ((edgeNum = triGeomShPtr->WhichEdge(edge)) > -1)
                            {
                                elementEdge = MemoryManager<ElementEdge>::AllocateSharedPtr();
                                elementEdge->m_Element = triGeomShPtr;
                                elementEdge->m_EdgeIndx = edgeNum;
                                returnval->push_back(elementEdge);
                            }
                        }
                        else if (quadGeomShPtr)
                        {
                            if ((edgeNum = quadGeomShPtr->WhichEdge(edge)) > -1)
                            {
                                elementEdge = MemoryManager<ElementEdge>::AllocateSharedPtr();
                                elementEdge->m_Element = quadGeomShPtr;
                                elementEdge->m_EdgeIndx = edgeNum;
                                returnval->push_back(elementEdge);
                            }
                        }
                    }
                }         
            }

            //for(triIter = m_trigeoms.begin(); triIter != m_trigeoms.end(); ++triIter)
            //{
            //    int edgeNum;
            //    if ((edgeNum = (*triIter)->WhichEdge(edge)) > -1)
            //    {
            //        elementEdge = MemoryManager<ElementEdge>::AllocateSharedPtr();
            //        elementEdge->m_Element = *triIter;
            //        elementEdge->m_EdgeIndx = edgeNum;
            //        returnval->push_back(elementEdge);
            //    }
            //}

            //QuadGeomVector::iterator quadIter;

            //for(quadIter = m_quadgeoms.begin(); quadIter != m_quadgeoms.end(); ++quadIter)
            //{
            //    int edgeNum;
            //    if ((edgeNum = (*quadIter)->WhichEdge(edge)) > -1)
            //    {
            //        elementEdge = MemoryManager<ElementEdge>::AllocateSharedPtr();
            //        elementEdge->m_Element = *quadIter;
            //        elementEdge->m_EdgeIndx = edgeNum;
            //        returnval->push_back(elementEdge);
            //    }
            //}

            return returnval;
        }

        LibUtilities::BasisKey MeshGraph2D::GetEdgeBasisKey(SegGeomSharedPtr edge)
        {
            ElementEdgeVectorSharedPtr elements = GetElementsFromEdge(edge);
            // Perhaps, a check should be done here to ensure that in case 
            // elements->size!=1, all elements to which the edge belongs have the same type
            // and order of expansion such that no confusion can arise.
            ExpansionShPtr expansion = GetExpansion((*elements)[0]->m_Element);
            

            int edge_id = (*elements)[0]->m_EdgeIndx;

            if((*elements)[0]->m_Element->GetGeomShapeType() == eTriangle)
            {
                edge_id = (edge_id)? 1:0;
            }
            else
            {
                edge_id = edge_id%2;
            }

            int nummodes  = expansion->m_BasisKeyVector[edge_id].GetNumModes();
            int numpoints = expansion->m_BasisKeyVector[edge_id].GetNumPoints();
            
	        if((*elements)[0]->m_Element->GetGeomShapeType() == eTriangle)
	        {
		        // Use edge 0 to define basis of order relevant to edge 
                switch(expansion->m_BasisKeyVector[edge_id].GetBasisType())
		        {
		            case LibUtilities::eModified_B:
                    {
		                switch(expansion->m_BasisKeyVector[edge_id].GetPointsType())
                        {
                            case LibUtilities::eGaussRadauMAlpha1Beta0:
			                {                            
			                    const LibUtilities::PointsKey pkey(numpoints+1,LibUtilities::eGaussLobattoLegendre);
			                    return LibUtilities::BasisKey(expansion->m_BasisKeyVector[0].GetBasisType(),nummodes,pkey);
			                }
			                break;

			                default:
			                ASSERTL0(false,"Unexpected points distribution");

                            // It doesn't matter what we return here since the ASSERT will stop execution.
                            // Just return something to prevent warnings messages.
                            const LibUtilities::PointsKey pkey(numpoints+1,LibUtilities::eGaussLobattoLegendre);
                            return LibUtilities::BasisKey(expansion->m_BasisKeyVector[0].GetBasisType(),nummodes,pkey);
			                break;
			            }
		            }
		            break;

		            case LibUtilities::eModified_A:
                    {
		                switch(expansion->m_BasisKeyVector[edge_id].GetPointsType())
                        {
                            case LibUtilities::eGaussLobattoLegendre:
			                {                            
			                    const LibUtilities::PointsKey pkey(numpoints,LibUtilities::eGaussLobattoLegendre);
			                    return LibUtilities::BasisKey(expansion->m_BasisKeyVector[0].GetBasisType(),nummodes,pkey);
			                }
			                break;

			                default:
			                ASSERTL0(false,"Unexpected points distribution");
                            // It doesn't matter what we return here since the ASSERT will stop execution.
                            // Just return something to prevent warnings messages.
                            const LibUtilities::PointsKey pkey(numpoints+1,LibUtilities::eGaussLobattoLegendre);
                            return LibUtilities::BasisKey(expansion->m_BasisKeyVector[0].GetBasisType(),nummodes,pkey);
			                break;
			            }
		            }
		            break;

		            default:
		            ASSERTL0(false,"Unexpected basis distribution");
                    // It doesn't matter what we return here since the ASSERT will stop execution.
                    // Just return something to prevent warnings messages.
                    const LibUtilities::PointsKey pkey(numpoints+1,LibUtilities::eGaussLobattoLegendre);
                    return LibUtilities::BasisKey(expansion->m_BasisKeyVector[0].GetBasisType(),nummodes,pkey);
		        }
	        }
	        else
	        {
		        // Quadrilateral
		        const LibUtilities::PointsKey pkey(numpoints,expansion->m_BasisKeyVector[edge_id].GetPointsType());
		        return LibUtilities::BasisKey(expansion->m_BasisKeyVector[edge_id].GetBasisType(),nummodes,pkey);
	        }
        }
    }; //end of namespace
}; //end of namespace

//
// $Log: MeshGraph2D.cpp,v $
// Revision 1.40  2009/10/22 17:41:47  cbiotto
// Update for variable order expansion
//
// Revision 1.39  2009/07/02 13:32:24  sehunchun
// *** empty log message ***
//
// Revision 1.38  2009/07/02 11:39:44  sehunchun
// Modification for 2D gemoetry embedded in 3D
//
{'fraction_non_alphanumeric': 0.10424592391304348, 'fraction_numerical': 0.005995244565217391, 'mean_word_length': 1.308968275753892, 'pattern_counts': {'":': 0, '<': 59, '<?xml version=': 0, '>': 262, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1067532', 'n_tokens_mistral': 15197, 'n_tokens_neox': 14033, 'n_words': 3255}