MeshGraph.cpp 171 KB
Newer Older
Mike Kirby's avatar
Mike Kirby committed
1 2
////////////////////////////////////////////////////////////////////////////////
//
Chris Cantwell's avatar
Chris Cantwell committed
3
//  File: MeshGraph.cpp
Mike Kirby's avatar
Mike Kirby committed
4 5 6 7 8 9
//
//  For more information, please see: http://www.nektar.info/
//
//  The MIT License
//
//  Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10
//  Department of Aeronautics, Imperial College London (UK), and Scientific
Mike Kirby's avatar
Mike Kirby committed
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
//  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.
//
32
//  Description:
Mike Kirby's avatar
Mike Kirby committed
33 34 35
//
////////////////////////////////////////////////////////////////////////////////

Joe Frazier's avatar
Joe Frazier committed
36

Mike Kirby's avatar
Mike Kirby committed
37
#include <SpatialDomains/MeshGraph.h>
38
#include <LibUtilities/BasicUtils/CompressData.h>
39
#include <LibUtilities/BasicUtils/ParseUtils.hpp>
40
#include <LibUtilities/BasicUtils/Equation.h>
41
#include <LibUtilities/BasicUtils/FieldIOXml.h>
42 43 44 45
#include <StdRegions/StdTriExp.h>
#include <StdRegions/StdTetExp.h>
#include <StdRegions/StdPyrExp.h>
#include <StdRegions/StdPrismExp.h>
Mike Kirby's avatar
Mike Kirby committed
46 47

// Use the stl version, primarily for string.
Blake Nelson's avatar
Blake Nelson committed
48
#ifndef TIXML_USE_STL
Mike Kirby's avatar
Mike Kirby committed
49
#define TIXML_USE_STL
Blake Nelson's avatar
Blake Nelson committed
50 51
#endif

52
#include <tinyxml.h>
53
#include <cstring>
Sophia Han's avatar
Sophia Han committed
54
#include <sstream>
55
#include <iomanip>
Mike Kirby's avatar
Mike Kirby committed
56

57 58
#include <SpatialDomains/MeshGraph1D.h>
#include <SpatialDomains/MeshGraph2D.h>
Michael DeLisi's avatar
Michael DeLisi committed
59
#include <SpatialDomains/MeshGraph3D.h>
60

61 62 63 64 65
// 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>
66
#include <boost/iostreams/filter/zlib.hpp>
67
#include <boost/iostreams/filtering_stream.hpp>
68
#include <boost/make_shared.hpp>
69

70 71
using namespace std;

Mike Kirby's avatar
Mike Kirby committed
72 73
namespace Nektar
{
74 75
    namespace SpatialDomains
    {
Chris Cantwell's avatar
Chris Cantwell committed
76 77 78
        /**
         *
         */
79 80
        MeshGraph::MeshGraph():
            m_meshDimension(3),
81 82
            m_spaceDimension(3),
            m_domainRange(NullDomainRangeShPtr)
83 84
        {
        }
85

Chris Cantwell's avatar
Chris Cantwell committed
86 87 88 89 90 91 92

        /**
         *
         */
        MeshGraph::MeshGraph(
                unsigned int meshDimension,
                unsigned int spaceDimension) :
93
            m_meshDimension(meshDimension),
94 95
            m_spaceDimension(spaceDimension),
            m_domainRange(NullDomainRangeShPtr)
96 97
        {
        }
98

Chris Cantwell's avatar
Chris Cantwell committed
99 100 101 102 103

        /**
         *
         */
        MeshGraph::MeshGraph(
104 105
                       const LibUtilities::SessionReaderSharedPtr &pSession,
                       const DomainRangeShPtr &rng):
106 107
            m_session(pSession),
            m_domainRange(rng)
108
        {
Chris Cantwell's avatar
Chris Cantwell committed
109
        }
110

Chris Cantwell's avatar
Chris Cantwell committed
111

112

Chris Cantwell's avatar
Chris Cantwell committed
113 114 115 116 117
        /**
         *
         */
        MeshGraph::~MeshGraph()
        {
118 119
        }

Chris Cantwell's avatar
Chris Cantwell committed
120 121 122 123

        /**
         *
         */
124
        boost::shared_ptr<MeshGraph> MeshGraph::Read(
Chris Cantwell's avatar
Chris Cantwell committed
125
                      const LibUtilities::SessionReaderSharedPtr &pSession,
126
                      DomainRangeShPtr &rng)
127 128 129 130 131 132 133 134 135 136 137 138 139
        {
            boost::shared_ptr<MeshGraph> returnval;

            // read the geometry tag to get the dimension

            TiXmlElement* geometry_tag = pSession->GetElement("NEKTAR/GEOMETRY");
            TiXmlAttribute *attr = geometry_tag->FirstAttribute();
            int meshDim = 0;
            while (attr)
            {
                std::string attrName(attr->Name());
                if (attrName == "DIM")
                {
140
                    int err = attr->QueryIntValue(&meshDim);
141
                    ASSERTL0(err==TIXML_SUCCESS, "Unable to read mesh dimension.");
142 143 144 145 146 147
                    break;
                }
                else
                {
                    std::string errstr("Unknown attribute: ");
                    errstr += attrName;
148
                    ASSERTL0(false, errstr.c_str());
149 150 151 152 153 154 155 156 157 158 159
                }

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

            // instantiate the dimension-specific meshgraph classes

            switch(meshDim)
            {
            case 1:
160
                returnval = MemoryManager<MeshGraph1D>::AllocateSharedPtr(pSession,rng);
161 162 163
                break;

            case 2:
164
                returnval = MemoryManager<MeshGraph2D>::AllocateSharedPtr(pSession,rng);
165 166 167
                break;

            case 3:
168
                returnval = MemoryManager<MeshGraph3D>::AllocateSharedPtr(pSession,rng);
169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184
                break;

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

            return returnval;
        }



        /*  ====  OUTDATED ROUTINE, PLEASE NOT USE  ==== */
Chris Cantwell's avatar
Chris Cantwell committed
185 186 187
        boost::shared_ptr<MeshGraph> MeshGraph::Read(
                const std::string& infilename,
                bool pReadExpansions)
188 189
        {
            boost::shared_ptr<MeshGraph> returnval;
190

191
            MeshGraph mesh;
192

193 194
            mesh.ReadGeometry(infilename);
            int meshDim = mesh.GetMeshDimension();
195

196 197 198 199 200
            switch(meshDim)
            {
            case 1:
                returnval = MemoryManager<MeshGraph1D>::AllocateSharedPtr();
                break;
201

202 203 204
            case 2:
                returnval = MemoryManager<MeshGraph2D>::AllocateSharedPtr();
                break;
205

206 207 208
            case 3:
                returnval = MemoryManager<MeshGraph3D>::AllocateSharedPtr();
                break;
209

210 211 212 213 214 215 216
            default:
                std::string err = "Invalid mesh dimension: ";
                std::stringstream strstrm;
                strstrm << meshDim;
                err += strstrm.str();
                NEKERROR(ErrorUtil::efatal, err.c_str());
            }
217

218 219 220 221
            if (returnval)
            {
                returnval->ReadGeometry(infilename);
                returnval->ReadGeometryInfo(infilename);
222 223 224 225
                if (pReadExpansions)
                {
                    returnval->ReadExpansions(infilename);
                }
226 227 228
            }
            return returnval;
        }
229

Chris Cantwell's avatar
Chris Cantwell committed
230

231

Chris Cantwell's avatar
Chris Cantwell committed
232 233 234 235
        /**
         *
         */
        void MeshGraph::ReadGeometry(const std::string& infilename)
236
        {
Chris Cantwell's avatar
Chris Cantwell committed
237 238
            TiXmlDocument doc(infilename);
            bool loadOkay = doc.LoadFile();
239

Chris Cantwell's avatar
Chris Cantwell committed
240 241 242 243 244
            std::stringstream errstr;
            errstr << "Unable to load file: " << infilename << " (";
            errstr << doc.ErrorDesc() << ", line " << doc.ErrorRow()
                                 << ", column " << doc.ErrorCol() << ")";
            ASSERTL0(loadOkay, errstr.str());
245

Chris Cantwell's avatar
Chris Cantwell committed
246 247
            ReadGeometry(doc);
        }
248 249


Chris Cantwell's avatar
Chris Cantwell committed
250 251 252 253 254 255 256 257
        /**
         *
         */
        void MeshGraph::ReadGeometry(TiXmlDocument &doc)
        {
            TiXmlHandle docHandle(&doc);
            TiXmlElement* mesh = NULL;
            TiXmlElement* master = NULL;    // Master tag within which all data is contained.
258

Chris Cantwell's avatar
Chris Cantwell committed
259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280
            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")
281
                {
Chris Cantwell's avatar
Chris Cantwell committed
282 283
                    err = attr->QueryIntValue(&m_meshDimension);
                    ASSERTL1(err==TIXML_SUCCESS, "Unable to read mesh dimension.");
284
                }
Chris Cantwell's avatar
Chris Cantwell committed
285
                else if (attrName == "SPACE")
286
                {
Chris Cantwell's avatar
Chris Cantwell committed
287 288 289 290 291 292 293 294 295 296 297 298 299 300 301
                    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());
                }
302

Chris Cantwell's avatar
Chris Cantwell committed
303 304 305
                // Get the next attribute.
                attr = attr->Next();
            }
306

Chris Cantwell's avatar
Chris Cantwell committed
307
            ASSERTL1(m_meshDimension<=m_spaceDimension, "Mesh dimension greater than space dimension");
308

Chris Cantwell's avatar
Chris Cantwell committed
309 310 311
            // Now read the vertices
            TiXmlElement* element = mesh->FirstChildElement("VERTEX");
            ASSERTL0(element, "Unable to find mesh VERTEX tag in file.");
312

Chris Cantwell's avatar
Chris Cantwell committed
313
            NekDouble xscale,yscale,zscale;
314

Chris Cantwell's avatar
Chris Cantwell committed
315 316 317 318 319 320 321 322 323 324 325 326 327 328 329
            // check to see if any scaling parameters are in
            // attributes and determine these values
            LibUtilities::AnalyticExpressionEvaluator expEvaluator;
            //LibUtilities::ExpressionEvaluator expEvaluator;
            const char *xscal =  element->Attribute("XSCALE");
            if(!xscal)
            {
                xscale = 1.0;
            }
            else
            {
                std::string xscalstr = xscal;
                int expr_id = expEvaluator.DefineFunction("",xscalstr);
                xscale = expEvaluator.Evaluate(expr_id);
            }
330

Chris Cantwell's avatar
Chris Cantwell committed
331 332 333 334 335 336 337 338 339 340 341
            const char *yscal =  element->Attribute("YSCALE");
            if(!yscal)
            {
                yscale = 1.0;
            }
            else
            {
                std::string yscalstr = yscal;
                int expr_id = expEvaluator.DefineFunction("",yscalstr);
                yscale = expEvaluator.Evaluate(expr_id);
            }
342

Chris Cantwell's avatar
Chris Cantwell committed
343 344 345 346 347 348 349 350 351 352 353
            const char *zscal = element->Attribute("ZSCALE");
            if(!zscal)
            {
                zscale = 1.0;
            }
            else
            {
                std::string zscalstr = zscal;
                int expr_id = expEvaluator.DefineFunction("",zscalstr);
                zscale = expEvaluator.Evaluate(expr_id);
            }
354

355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397

            NekDouble xmove,ymove,zmove;

            // check to see if any moving parameters are in
            // attributes and determine these values

            //LibUtilities::ExpressionEvaluator expEvaluator;
            const char *xmov =  element->Attribute("XMOVE");
            if(!xmov)
            {
                xmove = 0.0;
            }
            else
            {
                std::string xmovstr = xmov;
                int expr_id = expEvaluator.DefineFunction("",xmovstr);
                xmove = expEvaluator.Evaluate(expr_id);
            }

            const char *ymov =  element->Attribute("YMOVE");
            if(!ymov)
            {
                ymove = 0.0;
            }
            else
            {
                std::string ymovstr = ymov;
                int expr_id = expEvaluator.DefineFunction("",ymovstr);
                ymove = expEvaluator.Evaluate(expr_id);
            }

            const char *zmov = element->Attribute("ZMOVE");
            if(!zmov)
            {
                zmove = 0.0;
            }
            else
            {
                std::string zmovstr = zmov;
                int expr_id = expEvaluator.DefineFunction("",zmovstr);
                zmove = expEvaluator.Evaluate(expr_id);
            }

398 399 400 401
            string IsCompressed;
            element->QueryStringAttribute("COMPRESSED",&IsCompressed); 

            if(IsCompressed.size()) 
Chris Cantwell's avatar
Chris Cantwell committed
402
            {
403 404
                if(boost::iequals(IsCompressed,
                            LibUtilities::CompressData::GetCompressString()))
405
                {
406 407
                    // Extract the vertex body
                    TiXmlNode* vertexChild = element->FirstChild();
408 409 410 411
                    ASSERTL0(vertexChild,
                             "Unable to extract the data from the compressed "
                             "vertex tag.");

412 413 414 415 416
                    std::string vertexStr;
                    if (vertexChild->Type() == TiXmlNode::TINYXML_TEXT)
                    {
                        vertexStr += vertexChild->ToText()->ValueStr();
                    }
417

418
                    std::vector<LibUtilities::MeshVertex> vertData;
419 420 421
                    LibUtilities::CompressData::ZlibDecodeFromBase64Str(
                                                        vertexStr,vertData);

422 423 424 425 426 427 428 429
                    int indx;
                    NekDouble xval, yval, zval;
                    for(int i = 0; i < vertData.size(); ++i)
                    {
                        indx = vertData[i].id;
                        xval = vertData[i].x;
                        yval = vertData[i].y;
                        zval = vertData[i].z;
430

431 432 433
                        xval = xval*xscale + xmove;
                        yval = yval*yscale + ymove;
                        zval = zval*zscale + zmove;
434 435 436 437 438

                        PointGeomSharedPtr vert(
                            MemoryManager<PointGeom>::AllocateSharedPtr(
                                m_spaceDimension, indx, xval, yval, zval));

439 440 441
                        vert->SetGlobalID(indx);
                        m_vertSet[indx] = vert;
                    }
442
                }
443 444
                else
                {
445 446 447
                    ASSERTL0(false,"Compressed formats do not match. Expected :"
                             + LibUtilities::CompressData::GetCompressString()
                             + " but got " + std::string(IsCompressed));
448
                }
449 450 451 452
            }
            else
            {
                TiXmlElement *vertex = element->FirstChildElement("V");
453
                
454 455 456 457
                int indx;
                int nextVertexNumber = -1;
                
                while (vertex)
Chris Cantwell's avatar
Chris Cantwell committed
458
                {
459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474
                    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.");

                    // Now read body of vertex
                    std::string vertexBodyStr;
                    
                    TiXmlNode *vertexBody = vertex->FirstChild();
                    
                    while (vertexBody)
Chris Cantwell's avatar
Chris Cantwell committed
475
                    {
476 477
                        // Accumulate all non-comment body data.
                        if (vertexBody->Type() == TiXmlNode::TINYXML_TEXT)
478
                        {
479 480
                            vertexBodyStr += vertexBody->ToText()->Value();
                            vertexBodyStr += " ";
481
                        }
482 483
                        
                        vertexBody = vertexBody->NextSibling();
484
                    }
485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518
                    
                    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 = xval*xscale + xmove;
                            yval = yval*yscale + ymove;
                            zval = zval*zscale + zmove;
                            
                            // Need to check it here because we may not be
                            // good after the read indicating that there
                            // was nothing to read.
                            if (!vertexDataStrm.fail())
                            {
                                PointGeomSharedPtr vert(MemoryManager<PointGeom>::AllocateSharedPtr(m_spaceDimension, indx, xval, yval, zval));
                                vert->SetGlobalID(indx);
                                m_vertSet[indx] = vert;
                            }
                        }
                    }
                    catch(...)
                    {
                        ASSERTL0(false, "Unable to read VERTEX data.");
                    }
                    
                    vertex = vertex->NextSiblingElement("V");
519 520
                }
            }
521
        }
522 523


Chris Cantwell's avatar
Chris Cantwell committed
524 525 526 527 528 529 530
        /**
         * Read the geometry-related information from the given file. This
         * information is located within the XML tree under
         * <NEKTAR><GEOMETRY><GEOMINFO>.
         * @param   infilename      Filename of XML file.
         */
        void MeshGraph::ReadGeometryInfo(const std::string &infilename)
531
        {
Chris Cantwell's avatar
Chris Cantwell committed
532 533
            TiXmlDocument doc(infilename);
            bool loadOkay = doc.LoadFile();
534

Chris Cantwell's avatar
Chris Cantwell committed
535 536 537 538 539
            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());
540

Chris Cantwell's avatar
Chris Cantwell committed
541 542
            ReadGeometryInfo(doc);
        }
543 544


Chris Cantwell's avatar
Chris Cantwell committed
545 546 547 548 549 550 551 552 553 554
        /**
         * Read the geometry-related information from the given XML document.
         * This information is located within the XML tree under
         * <NEKTAR><GEOMETRY><GEOMINFO>.
         * @param   doc             XML document.
         */
        void MeshGraph::ReadGeometryInfo(TiXmlDocument &doc)
        {
            TiXmlElement *master = doc.FirstChildElement("NEKTAR");
            ASSERTL0(master, "Unable to find NEKTAR tag in file.");
555

Chris Cantwell's avatar
Chris Cantwell committed
556 557 558
            // Find the Expansions tag
            TiXmlElement *geomTag = master->FirstChildElement("GEOMETRY");
            ASSERTL0(geomTag, "Unable to find GEOMETRY tag in file.");
559

Chris Cantwell's avatar
Chris Cantwell committed
560 561 562
            // See if we have GEOMINFO. If there is none, it's fine.
            TiXmlElement *geomInfoTag = geomTag->FirstChildElement("GEOMINFO");
            if (!geomInfoTag) return;
563

Chris Cantwell's avatar
Chris Cantwell committed
564
            TiXmlElement *infoItem = geomInfoTag->FirstChildElement("I");
565

Chris Cantwell's avatar
Chris Cantwell committed
566 567 568 569 570 571 572
            // Multiple nodes will only occur if there is a comment in between
            // definitions.
            while (infoItem)
            {
                std::string geomProperty = infoItem->Attribute("PROPERTY");
                std::string geomValue    = infoItem->Attribute("VALUE");
                GeomInfoMap::iterator x  = m_geomInfo.find(geomProperty);
573

Chris Cantwell's avatar
Chris Cantwell committed
574 575 576 577 578 579
                ASSERTL0(x == m_geomInfo.end(),
                        "Property " + geomProperty + " already specified.");
                m_geomInfo[geomProperty] = geomValue;
                infoItem = infoItem->NextSiblingElement("I");
            }
        }
580 581


Chris Cantwell's avatar
Chris Cantwell committed
582 583 584 585 586 587 588
        /**
         *
         */
        void MeshGraph::ReadExpansions(const std::string& infilename)
        {
            TiXmlDocument doc(infilename);
            bool loadOkay = doc.LoadFile();
589

Chris Cantwell's avatar
Chris Cantwell committed
590 591 592 593 594
            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());
595

Chris Cantwell's avatar
Chris Cantwell committed
596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643
            ReadExpansions(doc);
        }


        /**
         *
         */
        void MeshGraph::ReadExpansions(TiXmlDocument &doc)
        {
            TiXmlElement *master = doc.FirstChildElement("NEKTAR");
            ASSERTL0(master, "Unable to find NEKTAR tag in file.");

            // Find the Expansions tag
            TiXmlElement *expansionTypes = master->FirstChildElement("EXPANSIONS");
            ASSERTL0(expansionTypes, "Unable to find EXPANSIONS tag in file.");

            if(expansionTypes)
            {
                // Find the Expansion type
                TiXmlElement *expansion = expansionTypes->FirstChildElement();
                std::string   expType   = expansion->Value();

                if(expType == "E")
                {
                    int i;
                    ExpansionMapShPtr expansionMap;

                    /// Expansiontypes will contain composite,
                    /// nummodes, and expansiontype (eModified, or
                    /// eOrthogonal) Or a full list of data of
                    /// basistype, nummodes, pointstype, numpoints;

                    /// Expansiontypes may also contain a list of
                    /// fields that this expansion relates to. If this
                    /// does not exist the variable is only set to
                    /// "DefaultVar".

                    while (expansion)
                    {

                        const char *fStr = expansion->Attribute("FIELDS");
                        std::vector<std::string> fieldStrings;

                        if(fStr) // extract other fields.
                        {
                            std::string fieldStr = fStr;
                            bool  valid = ParseUtils::GenerateOrderedStringVector(fieldStr.c_str(),fieldStrings);
                            ASSERTL0(valid,"Unable to correctly parse the field string in ExpansionTypes.");
644
                        }
Chris Cantwell's avatar
Chris Cantwell committed
645 646 647 648 649

                        // check to see if m_expasionVectorShPtrMap has
                        // already been intiailised and if not intiailse
                        // vector.
                        if(m_expansionMapShPtrMap.count("DefaultVar") == 0) // no previous definitions
650
                        {
Chris Cantwell's avatar
Chris Cantwell committed
651
                            expansionMap = SetUpExpansionMap();
652

Chris Cantwell's avatar
Chris Cantwell committed
653
                            m_expansionMapShPtrMap["DefaultVar"] = expansionMap;
654

Chris Cantwell's avatar
Chris Cantwell committed
655 656 657
                            // make sure all fields in this search point
                            // to same expansion vector;
                            for(i = 0; i < fieldStrings.size(); ++i)
658
                            {
Chris Cantwell's avatar
Chris Cantwell committed
659
                                m_expansionMapShPtrMap[fieldStrings[i]] = expansionMap;
660 661
                            }
                        }
Chris Cantwell's avatar
Chris Cantwell committed
662
                        else // default variable is defined
663
                        {
664

Chris Cantwell's avatar
Chris Cantwell committed
665
                            if(fieldStrings.size()) // fields are defined
666
                            {
Chris Cantwell's avatar
Chris Cantwell committed
667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688
                                //see if field exists
                                if(m_expansionMapShPtrMap.count(fieldStrings[0]))
                                {
                                    expansionMap = m_expansionMapShPtrMap.find(fieldStrings[0])->second;
                                }
                                else
                                {
                                    expansionMap = SetUpExpansionMap();
                                    // make sure all fields in this search point
                                    // to same expansion vector;
                                    for(i = 0; i < fieldStrings.size(); ++i)
                                    {
                                        if(m_expansionMapShPtrMap.count(fieldStrings[i]) == 0)
                                        {
                                            m_expansionMapShPtrMap[fieldStrings[i]] = expansionMap;
                                        }
                                        else
                                        {
                                            ASSERTL0(false,"Expansion vector for this field is already  setup");
                                        }
                                    }
                                }
689
                            }
Chris Cantwell's avatar
Chris Cantwell committed
690
                            else // use default variable list
691
                            {
Chris Cantwell's avatar
Chris Cantwell committed
692
                                expansionMap = m_expansionMapShPtrMap.find("DefaultVar")->second;
693
                            }
694 695

                        }
696

Chris Cantwell's avatar
Chris Cantwell committed
697 698 699 700 701 702
                        /// Mandatory components...optional are to follow later.
                        std::string compositeStr = expansion->Attribute("COMPOSITE");
                        ASSERTL0(compositeStr.length() > 3, "COMPOSITE must be specified in expansion definition");
                        int beg = compositeStr.find_first_of("[");
                        int end = compositeStr.find_first_of("]");
                        std::string compositeListStr = compositeStr.substr(beg+1,end-beg-1);
703

Chris Cantwell's avatar
Chris Cantwell committed
704 705
                        CompositeMap compositeVector;
                        GetCompositeList(compositeListStr, compositeVector);
706

Chris Cantwell's avatar
Chris Cantwell committed
707 708 709
                        bool          useExpansionType = false;
                        ExpansionType expansion_type;
                        int           num_modes;
710

Chris Cantwell's avatar
Chris Cantwell committed
711 712
                        LibUtilities::BasisKeyVector basiskeyvec;
                        const char * tStr = expansion->Attribute("TYPE");
713

Chris Cantwell's avatar
Chris Cantwell committed
714 715 716 717 718 719
                        if(tStr) // use type string to define expansion
                        {
                            std::string typeStr = tStr;
                            const std::string* begStr = kExpansionTypeStr;
                            const std::string* endStr = kExpansionTypeStr+eExpansionTypeSize;
                            const std::string* expStr = std::find(begStr, endStr, typeStr);
720

Chris Cantwell's avatar
Chris Cantwell committed
721 722
                            ASSERTL0(expStr != endStr, "Invalid expansion type.");
                            expansion_type = (ExpansionType)(expStr - begStr);
723 724


725 726 727 728
                            /// \todo solvers break the pattern 'instantiate Session -> instantiate MeshGraph'
                            /// and parse command line arguments by themselves; one needs to unify command
                            /// line arguments handling.
                            /// Solvers tend to call MeshGraph::Read statically -> m_session
729
                            /// is not defined -> no info about command line arguments presented
730 731
                            /// ASSERTL0(m_session != 0, "One needs to instantiate SessionReader first");

732 733 734
                            const char *nStr = expansion->Attribute("NUMMODES");
                            ASSERTL0(nStr,"NUMMODES was not defined in EXPANSION section of input");
                            std::string nummodesStr = nStr;
735

736 737 738 739 740 741 742 743 744 745
                            // ASSERTL0(m_session,"Session should be defined to evaluate nummodes ");
                            if (m_session)
                            {
                                LibUtilities::Equation nummodesEqn(m_session, nummodesStr);
                                num_modes = (int) nummodesEqn.Evaluate();
                            }
                            else
                            {
                                num_modes = boost::lexical_cast<int>(nummodesStr);
                            }
746

Chris Cantwell's avatar
Chris Cantwell committed
747 748 749 750 751 752 753 754
                            useExpansionType = true;
                        }
                        else // assume expansion is defined individually
                        {
                            // Extract the attributes.
                            const char *bTypeStr = expansion->Attribute("BASISTYPE");
                            ASSERTL0(bTypeStr,"TYPE or BASISTYPE was not defined in EXPANSION section of input");
                            std::string basisTypeStr = bTypeStr;
755

Chris Cantwell's avatar
Chris Cantwell committed
756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776
                            // interpret the basis type string.
                            std::vector<std::string> basisStrings;
                            std::vector<LibUtilities::BasisType> basis;
                            bool valid = ParseUtils::GenerateOrderedStringVector(basisTypeStr.c_str(), basisStrings);
                            ASSERTL0(valid, "Unable to correctly parse the basis types.");
                            for (vector<std::string>::size_type i = 0; i < basisStrings.size(); i++)
                            {
                                valid = false;
                                for (unsigned int j = 0; j < LibUtilities::SIZE_BasisType; j++)
                                {
                                    if (LibUtilities::BasisTypeMap[j] == basisStrings[i])
                                    {
                                        basis.push_back((LibUtilities::BasisType) j);
                                        valid = true;
                                        break;
                                    }
                                }
                                ASSERTL0(valid, std::string("Unable to correctly parse the basis type: ").append(basisStrings[i]).c_str());
                            }
                            const char *nModesStr = expansion->Attribute("NUMMODES");
                            ASSERTL0(nModesStr,"NUMMODES was not defined in EXPANSION section of input");
777

Chris Cantwell's avatar
Chris Cantwell committed
778 779 780 781 782
                            std::string numModesStr = nModesStr;
                            std::vector<unsigned int> numModes;
                            valid = ParseUtils::GenerateOrderedVector(numModesStr.c_str(), numModes);
                            ASSERTL0(valid, "Unable to correctly parse the number of modes.");
                            ASSERTL0(numModes.size() == basis.size(),"information for num modes does not match the number of basis");
783

Chris Cantwell's avatar
Chris Cantwell committed
784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805
                            const char *pTypeStr =  expansion->Attribute("POINTSTYPE");
                            ASSERTL0(pTypeStr,"POINTSTYPE was not defined in EXPANSION section of input");
                            std::string pointsTypeStr = pTypeStr;
                            // interpret the points type string.
                            std::vector<std::string> pointsStrings;
                            std::vector<LibUtilities::PointsType> points;
                            valid = ParseUtils::GenerateOrderedStringVector(pointsTypeStr.c_str(), pointsStrings);
                            ASSERTL0(valid, "Unable to correctly parse the points types.");
                            for (vector<std::string>::size_type i = 0; i < pointsStrings.size(); i++)
                            {
                                valid = false;
                                for (unsigned int j = 0; j < LibUtilities::SIZE_PointsType; j++)
                                {
                                    if (LibUtilities::kPointsTypeStr[j] == pointsStrings[i])
                                    {
                                        points.push_back((LibUtilities::PointsType) j);
                                        valid = true;
                                        break;
                                    }
                                }
                                ASSERTL0(valid, std::string("Unable to correctly parse the points type: ").append(pointsStrings[i]).c_str());
                            }
806

Chris Cantwell's avatar
Chris Cantwell committed
807 808 809 810 811 812 813
                            const char *nPointsStr = expansion->Attribute("NUMPOINTS");
                            ASSERTL0(nPointsStr,"NUMPOINTS was not defined in EXPANSION section of input");
                            std::string numPointsStr = nPointsStr;
                            std::vector<unsigned int> numPoints;
                            valid = ParseUtils::GenerateOrderedVector(numPointsStr.c_str(), numPoints);
                            ASSERTL0(valid, "Unable to correctly parse the number of points.");
                            ASSERTL0(numPoints.size() == numPoints.size(),"information for num points does not match the number of basis");
814

Chris Cantwell's avatar
Chris Cantwell committed
815 816 817 818 819 820 821
                            for(int i = 0; i < basis.size(); ++i)
                            {
                                //Generate Basis key  using information
                                const LibUtilities::PointsKey pkey(numPoints[i],points[i]);
                                basiskeyvec.push_back(LibUtilities::BasisKey(basis[i],numModes[i],pkey));
                            }
                        }
822

Chris Cantwell's avatar
Chris Cantwell committed
823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845
                        // Now have composite and basiskeys.  Cycle through
                        // all composites for the geomShPtrs and set the modes
                        // and types for the elements contained in the element
                        // list.
                        CompositeMapIter compVecIter;
                        for (compVecIter = compositeVector.begin(); compVecIter != compositeVector.end(); ++compVecIter)
                        {
                            GeometryVectorIter geomVecIter;
                            for (geomVecIter = (compVecIter->second)->begin(); geomVecIter != (compVecIter->second)->end(); ++geomVecIter)
                            {
                                ExpansionMapIter x = expansionMap->find((*geomVecIter)->GetGlobalID());
                                ASSERTL0(x != expansionMap->end(), "Expansion not found!!");
                                if(useExpansionType)
                                {
                                    (x->second)->m_basisKeyVector = MeshGraph::DefineBasisKeyFromExpansionType(*geomVecIter,expansion_type,num_modes);
                                }
                                else
                                {
                                    ASSERTL0((*geomVecIter)->GetShapeDim() == basiskeyvec.size()," There is an incompatible expansion dimension with geometry dimension");
                                    (x->second)->m_basisKeyVector = basiskeyvec;
                                }
                            }
                        }
846

Chris Cantwell's avatar
Chris Cantwell committed
847 848 849 850 851 852 853
                        expansion = expansion->NextSiblingElement("E");
                    }
                }
                else if(expType == "H")
                {
                    int i;
                    ExpansionMapShPtr expansionMap;
854

Chris Cantwell's avatar
Chris Cantwell committed
855 856
                    while (expansion)
                    {
857

Chris Cantwell's avatar
Chris Cantwell committed
858 859
                        const char *fStr = expansion->Attribute("FIELDS");
                        std::vector<std::string> fieldStrings;
860

Chris Cantwell's avatar
Chris Cantwell committed
861 862 863 864 865 866
                        if(fStr) // extract other fields.
                        {
                            std::string fieldStr = fStr;
                            bool  valid = ParseUtils::GenerateOrderedStringVector(fieldStr.c_str(),fieldStrings);
                            ASSERTL0(valid,"Unable to correctly parse the field string in ExpansionTypes.");
                        }
867

Chris Cantwell's avatar
Chris Cantwell committed
868 869 870 871 872 873
                        // check to see if m_expasionVectorShPtrMap has
                        // already been intiailised and if not intiailse
                        // vector.
                        if(m_expansionMapShPtrMap.count("DefaultVar") == 0) // no previous definitions
                        {
                            expansionMap = SetUpExpansionMap();
874

Chris Cantwell's avatar
Chris Cantwell committed
875
                            m_expansionMapShPtrMap["DefaultVar"] = expansionMap;
876

Chris Cantwell's avatar
Chris Cantwell committed
877 878 879 880 881 882 883 884 885
                            // make sure all fields in this search point
                            // to same expansion vector;
                            for(i = 0; i < fieldStrings.size(); ++i)
                            {
                                m_expansionMapShPtrMap[fieldStrings[i]] = expansionMap;
                            }
                        }
                        else // default variable is defined
                        {
886

Chris Cantwell's avatar
Chris Cantwell committed
887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915
                            if(fieldStrings.size()) // fields are defined
                            {
                                //see if field exists
                                if(m_expansionMapShPtrMap.count(fieldStrings[0]))
                                {
                                    expansionMap = m_expansionMapShPtrMap.find(fieldStrings[0])->second;
                                }
                                else
                                {
                                    expansionMap = SetUpExpansionMap();
                                    // make sure all fields in this search point
                                    // to same expansion vector;
                                    for(i = 0; i < fieldStrings.size(); ++i)
                                    {
                                        if(m_expansionMapShPtrMap.count(fieldStrings[i]) == 0)
                                        {
                                            m_expansionMapShPtrMap[fieldStrings[i]] = expansionMap;
                                        }
                                        else
                                        {
                                            ASSERTL0(false,"Expansion vector for this field is already  setup");
                                        }
                                    }
                                }
                            }
                            else // use default variable list
                            {
                                expansionMap = m_expansionMapShPtrMap.find("DefaultVar")->second;
                            }
916

Chris Cantwell's avatar
Chris Cantwell committed
917
                        }
918

Chris Cantwell's avatar
Chris Cantwell committed
919 920 921 922 923 924
                        /// Mandatory components...optional are to follow later.
                        std::string compositeStr = expansion->Attribute("COMPOSITE");
                        ASSERTL0(compositeStr.length() > 3, "COMPOSITE must be specified in expansion definition");
                        int beg = compositeStr.find_first_of("[");
                        int end = compositeStr.find_first_of("]");
                        std::string compositeListStr = compositeStr.substr(beg+1,end-beg-1);
925

Chris Cantwell's avatar
Chris Cantwell committed
926 927
                        CompositeMap compositeVector;
                        GetCompositeList(compositeListStr, compositeVector);
928

929 930 931 932 933 934
                        ExpansionType expansion_type_x = eNoExpansionType;
                        ExpansionType expansion_type_y = eNoExpansionType;
                        ExpansionType expansion_type_z = eNoExpansionType;
                        int           num_modes_x = 0;
                        int           num_modes_y = 0;
                        int           num_modes_z = 0;
935

Chris Cantwell's avatar
Chris Cantwell committed
936
                        LibUtilities::BasisKeyVector basiskeyvec;
937

Chris Cantwell's avatar
Chris Cantwell committed
938
                        const char * tStr_x = expansion->Attribute("TYPE-X");
939

Chris Cantwell's avatar
Chris Cantwell committed
940
                        if(tStr_x) // use type string to define expansion
941
                        {
Chris Cantwell's avatar
Chris Cantwell committed
942 943 944 945
                            std::string typeStr = tStr_x;
                            const std::string* begStr = kExpansionTypeStr;
                            const std::string* endStr = kExpansionTypeStr+eExpansionTypeSize;
                            const std::string* expStr = std::find(begStr, endStr, typeStr);
946

Chris Cantwell's avatar
Chris Cantwell committed
947 948
                            ASSERTL0(expStr != endStr, "Invalid expansion type.");
                            expansion_type_x = (ExpansionType)(expStr - begStr);
949

Chris Cantwell's avatar
Chris Cantwell committed
950 951 952
                            const char *nStr = expansion->Attribute("NUMMODES-X");
                            ASSERTL0(nStr,"NUMMODES-X was not defined in EXPANSION section of input");
                            std::string nummodesStr = nStr;
953

954 955 956 957 958 959 960 961 962 963 964
                            // ASSERTL0(m_session,"Session should be defined to evaluate nummodes ");

                            if (m_session)
                            {
                                LibUtilities::Equation nummodesEqn(m_session, nummodesStr);
                                num_modes_x = (int) nummodesEqn.Evaluate();
                            }
                            else
                            {
                                num_modes_x = boost::lexical_cast<int>(nummodesStr);
                            }
965

Chris Cantwell's avatar
Chris Cantwell committed
966
                        }
967

Chris Cantwell's avatar
Chris Cantwell committed
968
                        const char * tStr_y = expansion->Attribute("TYPE-Y");
969

Chris Cantwell's avatar
Chris Cantwell committed
970 971 972 973 974 975 976 977 978 979 980 981 982 983
                        if(tStr_y) // use type string to define expansion
                        {
                            std::string typeStr = tStr_y;
                            const std::string* begStr = kExpansionTypeStr;
                            const std::string* endStr = kExpansionTypeStr+eExpansionTypeSize;
                            const std::string* expStr = std::find(begStr, endStr, typeStr);

                            ASSERTL0(expStr != endStr, "Invalid expansion type.");
                            expansion_type_y = (ExpansionType)(expStr - begStr);

                            const char *nStr = expansion->Attribute("NUMMODES-Y");
                            ASSERTL0(nStr,"NUMMODES-Y was not defined in EXPANSION section of input");
                            std::string nummodesStr = nStr;

984 985 986 987 988 989 990 991 992 993
                            // ASSERTL0(m_session,"Session should be defined to evaluate nummodes ");
                            if (m_session)
                            {
                                LibUtilities::Equation nummodesEqn(m_session, nummodesStr);
                                num_modes_y = (int) nummodesEqn.Evaluate();
                            }
                            else
                            {
                                num_modes_y = boost::lexical_cast<int>(nummodesStr);
                            }
Chris Cantwell's avatar
Chris Cantwell committed
994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012

                        }

                        const char * tStr_z = expansion->Attribute("TYPE-Z");

                        if(tStr_z) // use type string to define expansion
                        {
                            std::string typeStr = tStr_z;
                            const std::string* begStr = kExpansionTypeStr;
                            const std::string* endStr = kExpansionTypeStr+eExpansionTypeSize;
                            const std::string* expStr = std::find(begStr, endStr, typeStr);

                            ASSERTL0(expStr != endStr, "Invalid expansion type.");
                            expansion_type_z = (ExpansionType)(expStr - begStr);

                            const char *nStr = expansion->Attribute("NUMMODES-Z");
                            ASSERTL0(nStr,"NUMMODES-Z was not defined in EXPANSION section of input");
                            std::string nummodesStr = nStr;

1013 1014 1015 1016 1017 1018 1019 1020 1021 1022
                            // ASSERTL0(m_session,"Session should be defined to evaluate nummodes ");
                            if (m_session)
                            {
                                LibUtilities::Equation nummodesEqn(m_session, nummodesStr);
                                num_modes_z = (int) nummodesEqn.Evaluate();
                            }
                            else
                            {
                                num_modes_z = boost::lexical_cast<int>(nummodesStr);
                            }
Chris Cantwell's avatar
Chris Cantwell committed
1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051

                        }

                        CompositeMapIter compVecIter;
                        for (compVecIter = compositeVector.begin(); compVecIter != compositeVector.end(); ++compVecIter)
                        {
                            GeometryVectorIter geomVecIter;
                            for (geomVecIter = (compVecIter->second)->begin(); geomVecIter != (compVecIter->second)->end(); ++geomVecIter)
                            {
                                ExpansionMapIter expVecIter;
                                for (expVecIter = expansionMap->begin(); expVecIter != expansionMap->end(); ++expVecIter)
                                {

                                    (expVecIter->second)->m_basisKeyVector = DefineBasisKeyFromExpansionTypeHomo(*geomVecIter,
                                            expansion_type_x,
                                            expansion_type_y,
                                            expansion_type_z,
                                            num_modes_x,
                                            num_modes_y,
                                            num_modes_z);
                                }
                            }
                        }

                        expansion = expansion->NextSiblingElement("H");
                    }
                }
                else if(expType == "ELEMENTS")  // Reading a file with the expansion definition
                {
1052
                    std::vector<LibUtilities::FieldDefinitionsSharedPtr> fielddefs;
1053 1054

                    // This has to use the XML reader since we are treating the already parsed XML as a standard FLD file.
1055
                    boost::shared_ptr<LibUtilities::FieldIOXml> f = boost::make_shared<LibUtilities::FieldIOXml>(m_session->GetComm(), false);
1056
                    f->ImportFieldDefs(LibUtilities::XmlDataSource::create(doc), fielddefs, true);
Chris Cantwell's avatar
Chris Cantwell committed
1057 1058 1059
                    cout << "    Number of elements: " << fielddefs.size() << endl;
                    SetExpansions(fielddefs);
                }
1060 1061 1062 1063 1064 1065 1066 1067 1068 1069
                else if(expType == "F")
                {
                    ASSERTL0(expansion->Attribute("FILE"),
                                "Attribute FILE expected for type F expansion");
                    std::string filenameStr = expansion->Attribute("FILE");
                    ASSERTL0(!filenameStr.empty(),
                                "A filename must be specified for the FILE "
                                "attribute of expansion");

                    std::vector<LibUtilities::FieldDefinitionsSharedPtr> fielddefs;
1070 1071 1072 1073
                    LibUtilities::FieldIOSharedPtr f =
                            LibUtilities::FieldIO::CreateForFile(
                            m_session, filenameStr);
                    f->Import(filenameStr, fielddefs);
1074 1075
                    SetExpansions(fielddefs);
                }
Chris Cantwell's avatar
Chris Cantwell committed
1076 1077 1078 1079 1080
                else
                {
                    ASSERTL0(false,"Expansion type not defined");
                }
            }
1081
        }
1082

Chris Cantwell's avatar
Chris Cantwell committed
1083 1084 1085 1086 1087

        /**
         *
         */
        void MeshGraph::ReadDomain(TiXmlDocument &doc)
1088
        {
Chris Cantwell's avatar
Chris Cantwell committed
1089
            TiXmlHandle docHandle(&doc);
1090

Chris Cantwell's avatar
Chris Cantwell committed
1091 1092 1093 1094 1095 1096 1097 1098 1099 1100
            TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
            TiXmlElement* domain = NULL;

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

            /// Look for data in DOMAIN block.
            domain = mesh->FirstChildElement("DOMAIN");

            ASSERTL0(domain, "Unable to find DOMAIN tag in file.");

1101 1102 1103
            /// Elements are of the form: "<D ID = "N"> ... </D>".
            /// Read the ID field first.
            TiXmlElement *multidomains = domain->FirstChildElement("D");
1104

1105 1106 1107 1108 1109 1110 1111
            if(multidomains)
            {
                int nextDomainNumber = 0;
                while (multidomains)
                {
                    int indx;
                    int err = multidomains->QueryIntAttribute("ID", &indx);
Chris Cantwell's avatar
Chris Cantwell committed
1112
                    ASSERTL0(err == TIXML_SUCCESS,
1113
                             "Unable to read attribute ID in Domain.");
1114 1115


1116
                    TiXmlNode* elementChild = multidomains->FirstChild();
Dave Moxey's avatar
Dave Moxey committed
1117
                    while(elementChild && elementChild->Type() != TiXmlNode::TINYXML_TEXT)
1118 1119 1120
                    {
                        elementChild = elementChild->NextSibling();
                    }
Chris Cantwell's avatar
Chris Cantwell committed
1121

1122 1123
                    ASSERTL0(elementChild, "Unable to read DOMAIN body.");
                    std::string elementStr = elementChild->ToText()->ValueStr();
Chris Cantwell's avatar
Chris Cantwell committed
1124

1125
                    elementStr = elementStr.substr(elementStr.find_first_not_of(" "));
Chris Cantwell's avatar
Chris Cantwell committed
1126

1127 1128 1129
                    std::string::size_type indxBeg = elementStr.find_first_of('[') + 1;
                    std::string::size_type indxEnd = elementStr.find_last_of(']') - 1;
                    std::string indxStr = elementStr.substr(indxBeg, indxEnd - indxBeg + 1);
Chris Cantwell's avatar
Chris Cantwell committed
1130

1131
                    ASSERTL0(!indxStr.empty(), "Unable to read domain's composite index (index missing?).");
Chris Cantwell's avatar
Chris Cantwell committed
1132

1133 1134 1135 1136 1137
                    // Read the domain composites.
                    // Parse the composites into a list.
                    CompositeMap unrollDomain;
                    GetCompositeList(indxStr, unrollDomain);
                    m_domain.push_back(unrollDomain);
Chris Cantwell's avatar
Chris Cantwell committed
1138

1139 1140 1141 1142 1143
                    ASSERTL0(!m_domain[nextDomainNumber++].empty(), (std::string("Unable to obtain domain's referenced composite: ") + indxStr).c_str());

                    /// Keep looking
                    multidomains = multidomains->NextSiblingElement("D");
                }
1144

1145 1146 1147
            }
            else // previous definition of just one composite
            {
1148

1149 1150
                // find the non comment portion of the body.
                TiXmlNode* elementChild = domain->FirstChild();
Dave Moxey's avatar
Dave Moxey committed
1151
                while(elementChild && elementChild->Type() != TiXmlNode::TINYXML_TEXT)
1152 1153 1154
                {
                    elementChild = elementChild->NextSibling();
                }
Chris Cantwell's avatar
Chris Cantwell committed
1155

1156 1157
                ASSERTL0(elementChild, "Unable to read DOMAIN body.");
                std::string elementStr = elementChild->ToText()->ValueStr();
Chris Cantwell's avatar
Chris Cantwell committed
1158

1159
                elementStr = elementStr.substr(elementStr.find_first_not_of(" "));
Chris Cantwell's avatar
Chris Cantwell committed
1160

1161 1162 1163
                std::string::size_type indxBeg = elementStr.find_first_of('[') + 1;
                std::string::size_type indxEnd = elementStr.find_last_of(']') - 1;
                std::string indxStr = elementStr.substr(indxBeg, indxEnd - indxBeg + 1);
Chris Cantwell's avatar
Chris Cantwell committed
1164

1165
                ASSERTL0(!indxStr.empty(), "Unable to read domain's composite index (index missing?).");
Chris Cantwell's avatar
Chris Cantwell committed
1166

1167 1168 1169 1170 1171
                // Read the domain composites.
                // Parse the composites into a list.
                CompositeMap fullDomain;
                GetCompositeList(indxStr, fullDomain);
                m_domain.push_back(fullDomain);
Chris Cantwell's avatar
Chris Cantwell committed
1172

1173 1174
                ASSERTL0(!m_domain[0].empty(), (std::string("Unable to obtain domain's referenced composite: ") + indxStr).c_str());
            }
1175
        }
1176

Chris Cantwell's avatar
Chris Cantwell committed
1177 1178 1179 1180 1181

        /**
         *
         */
        void MeshGraph::ReadCurves(TiXmlDocument &doc)
1182
        {
Chris Cantwell's avatar
Chris Cantwell committed
1183 1184 1185 1186
            /// 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;
1187

Chris Cantwell's avatar
Chris Cantwell committed
1188 1189 1190 1191
            // check to see if any scaling parameters are in
            // attributes and determine these values
            TiXmlElement* element = mesh->FirstChildElement("VERTEX");
            ASSERTL0(element, "Unable to find mesh VERTEX tag in file.");
1192

Chris Cantwell's avatar
Chris Cantwell committed
1193
            NekDouble xscale,yscale,zscale;
1194

Chris Cantwell's avatar
Chris Cantwell committed
1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206
            LibUtilities::AnalyticExpressionEvaluator expEvaluator;
            const char *xscal =  element->Attribute("XSCALE");
            if(!xscal)
            {
                xscale = 1.0;
            }
            else
            {
                std::string xscalstr = xscal;
                int expr_id = expEvaluator.DefineFunction("",xscalstr);
                xscale = expEvaluator.Evaluate(expr_id);
            }
1207

Chris Cantwell's avatar
Chris Cantwell committed
1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218
            const char *yscal =  element->Attribute("YSCALE");
            if(!yscal)
            {
                yscale = 1.0;
            }
            else
            {
                std::string yscalstr = yscal;
                int expr_id = expEvaluator.DefineFunction("",yscalstr);
                yscale = expEvaluator.Evaluate(expr_id);
            }
1219

Chris Cantwell's avatar
Chris Cantwell committed
1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230
            const char *zscal = element->Attribute("ZSCALE");
            if(!zscal)
            {
                zscale = 1.0;
            }
            else
            {
                std::string zscalstr = zscal;
                int expr_id = expEvaluator.DefineFunction("",zscalstr);
                zscale = expEvaluator.Evaluate(expr_id);
            }
1231

1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273
            NekDouble xmove,ymove,zmove;

            // check to see if any moving parameters are in
            // attributes and determine these values

            //LibUtilities::ExpressionEvaluator expEvaluator;
            const char *xmov =  element->Attribute("XMOVE");
            if(!xmov)
            {
                xmove = 0.0;
            }
            else
            {
                std::string xmovstr = xmov;
                int expr_id = expEvaluator.DefineFunction("",xmovstr);
                xmove = expEvaluator.Evaluate(expr_id);
            }

            const char *ymov =  element->Attribute("YMOVE");
            if(!ymov)
            {
                ymove = 0.0;
            }
            else
            {
                std::string ymovstr = ymov;
                int expr_id = expEvaluator.DefineFunction("",ymovstr);
                ymove = expEvaluator.Evaluate(expr_id);
            }

            const char *zmov = element->Attribute("ZMOVE");
            if(!zmov)
            {
                zmove = 0.0;
            }
            else
            {
                std::string zmovstr = zmov;
                int expr_id = expEvaluator.DefineFunction("",zmovstr);
                zmove = expEvaluator.Evaluate(expr_id);
            }

Chris Cantwell's avatar
Chris Cantwell committed
1274
            int err;
1275

Chris Cantwell's avatar
Chris Cantwell committed
1276 1277
            /// Look for elements in CURVE block.
            field = mesh->FirstChildElement("CURVED");
1278

Chris Cantwell's avatar
Chris Cantwell committed
1279 1280 1281 1282 1283
            if(!field) //return if no curved entities
            {
                return;
            }

1284 1285 1286 1287
            string IsCompressed;
            field->QueryStringAttribute("COMPRESSED",&IsCompressed); 
            
            if(IsCompressed.size()) 
1288
            {
1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301 1302
                ASSERTL0(boost::iequals(IsCompressed,
                            LibUtilities::CompressData::GetCompressString()),
                         "Compressed formats do not match. Expected :"
                         + LibUtilities::CompressData::GetCompressString()
                         + " but got "
                         + boost::lexical_cast<std::string>(IsCompressed));

                std::vector<LibUtilities::MeshCurvedInfo> edginfo;
                std::vector<LibUtilities::MeshCurvedInfo> facinfo;
                LibUtilities::MeshCurvedPts cpts;

                // read edge, face info and curved poitns.
                TiXmlElement *x = field->FirstChildElement();
                while(x)
1303
                {
1304 1305 1306
                    const char *entitytype = x->Value();
                    // read in edge or face info
                    if(boost::iequals(entitytype,"E"))
1307
                    {
1308 1309 1310 1311 1312
                        // read in data
                        std::string elmtStr;
                        TiXmlNode* child = x->FirstChild();

                        if (child->Type() == TiXmlNode::TINYXML_TEXT)
1313
                        {
1314
                            elmtStr += child->ToText()->ValueStr();
1315
                        }
1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326

                        LibUtilities::CompressData::ZlibDecodeFromBase64Str(
                                                        elmtStr,edginfo);
                    }
                    else if(boost::iequals(entitytype,"F"))
                    {
                        // read in data
                        std::string elmtStr;
                        TiXmlNode* child = x->FirstChild();

                        if (child->Type() == TiXmlNode::TINYXML_TEXT)
1327
                        {
1328 1329 1330 1331 1332
                            elmtStr += child->ToText()->ValueStr();
                        }

                        LibUtilities::CompressData::ZlibDecodeFromBase64Str(
                                                        elmtStr,facinfo);
1333
                    }
1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351
                    else if(boost::iequals(entitytype,"DATAPOINTS"))
                    {
                        NekInt id;
                        ASSERTL0(x->Attribute("ID", &id),
                                 "Failed to get ID from PTS section");
                        cpts.id = id;

                        // read in data
                        std::string elmtStr;

                        TiXmlElement* DataIdx =
                            x->FirstChildElement("INDEX");
                        ASSERTL0(DataIdx,
                                 "Cannot read data index tag in compressed "
                                 "curved section");

                        TiXmlNode* child = DataIdx->FirstChild();
                        if (child->Type() == TiXmlNode::TINYXML_TEXT)
1352
                        {
1353
                            elmtStr = child->ToText()->ValueStr();
1354
                        }