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

                        }

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

                        if(tStr_z) // use type string to define expansion
                        {