Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MeshGraph.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: MeshGraph.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description:
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 
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 #include <StdRegions/StdTriExp.h>
43 #include <StdRegions/StdTetExp.h>
44 #include <StdRegions/StdPyrExp.h>
45 #include <StdRegions/StdPrismExp.h>
46 
47 // Use the stl version, primarily for string.
48 #ifndef TIXML_USE_STL
49 #define TIXML_USE_STL
50 #endif
51 
52 #include <tinyxml.h>
53 #include <cstring>
54 #include <sstream>
55 #include <iomanip>
56 
57 #include <SpatialDomains/MeshGraph1D.h>
58 #include <SpatialDomains/MeshGraph2D.h>
59 #include <SpatialDomains/MeshGraph3D.h>
60 
61 // These are required for the Write(...) and Import(...) functions.
62 #include <boost/archive/iterators/base64_from_binary.hpp>
63 #include <boost/archive/iterators/binary_from_base64.hpp>
64 #include <boost/archive/iterators/transform_width.hpp>
65 #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 using namespace std;
71 
72 namespace Nektar
73 {
74  namespace SpatialDomains
75  {
76  /**
77  *
78  */
79  MeshGraph::MeshGraph():
80  m_meshDimension(3),
81  m_spaceDimension(3),
82  m_domainRange(NullDomainRangeShPtr)
83  {
84  }
85 
86 
87  /**
88  *
89  */
90  MeshGraph::MeshGraph(
91  unsigned int meshDimension,
92  unsigned int spaceDimension) :
93  m_meshDimension(meshDimension),
94  m_spaceDimension(spaceDimension),
95  m_domainRange(NullDomainRangeShPtr)
96  {
97  }
98 
99 
100  /**
101  *
102  */
103  MeshGraph::MeshGraph(
104  const LibUtilities::SessionReaderSharedPtr &pSession,
105  const DomainRangeShPtr &rng):
106  m_session(pSession),
107  m_domainRange(rng)
108  {
109  }
110 
111 
112 
113  /**
114  *
115  */
116  MeshGraph::~MeshGraph()
117  {
118  }
119 
120 
121  /**
122  *
123  */
124  boost::shared_ptr<MeshGraph> MeshGraph::Read(
125  const LibUtilities::SessionReaderSharedPtr &pSession,
126  DomainRangeShPtr &rng)
127  {
128  boost::shared_ptr<MeshGraph> returnval;
129 
130  // read the geometry tag to get the dimension
131 
132  TiXmlElement* geometry_tag = pSession->GetElement("NEKTAR/GEOMETRY");
133  TiXmlAttribute *attr = geometry_tag->FirstAttribute();
134  int meshDim = 0;
135  while (attr)
136  {
137  std::string attrName(attr->Name());
138  if (attrName == "DIM")
139  {
140  int err = attr->QueryIntValue(&meshDim);
141  ASSERTL0(err==TIXML_SUCCESS, "Unable to read mesh dimension.");
142  break;
143  }
144  else
145  {
146  std::string errstr("Unknown attribute: ");
147  errstr += attrName;
148  ASSERTL0(false, errstr.c_str());
149  }
150 
151  // Get the next attribute.
152  attr = attr->Next();
153  }
154 
155  // instantiate the dimension-specific meshgraph classes
156 
157  switch(meshDim)
158  {
159  case 1:
160  returnval = MemoryManager<MeshGraph1D>::AllocateSharedPtr(pSession,rng);
161  break;
162 
163  case 2:
164  returnval = MemoryManager<MeshGraph2D>::AllocateSharedPtr(pSession,rng);
165  break;
166 
167  case 3:
168  returnval = MemoryManager<MeshGraph3D>::AllocateSharedPtr(pSession,rng);
169  break;
170 
171  default:
172  std::string err = "Invalid mesh dimension: ";
173  std::stringstream strstrm;
174  strstrm << meshDim;
175  err += strstrm.str();
176  NEKERROR(ErrorUtil::efatal, err.c_str());
177  }
178 
179  return returnval;
180  }
181 
182 
183 
184  /* ==== OUTDATED ROUTINE, PLEASE NOT USE ==== */
185  boost::shared_ptr<MeshGraph> MeshGraph::Read(
186  const std::string& infilename,
187  bool pReadExpansions)
188  {
189  boost::shared_ptr<MeshGraph> returnval;
190 
191  MeshGraph mesh;
192 
193  mesh.ReadGeometry(infilename);
194  int meshDim = mesh.GetMeshDimension();
195 
196  switch(meshDim)
197  {
198  case 1:
199  returnval = MemoryManager<MeshGraph1D>::AllocateSharedPtr();
200  break;
201 
202  case 2:
203  returnval = MemoryManager<MeshGraph2D>::AllocateSharedPtr();
204  break;
205 
206  case 3:
207  returnval = MemoryManager<MeshGraph3D>::AllocateSharedPtr();
208  break;
209 
210  default:
211  std::string err = "Invalid mesh dimension: ";
212  std::stringstream strstrm;
213  strstrm << meshDim;
214  err += strstrm.str();
215  NEKERROR(ErrorUtil::efatal, err.c_str());
216  }
217 
218  if (returnval)
219  {
220  returnval->ReadGeometry(infilename);
221  returnval->ReadGeometryInfo(infilename);
222  if (pReadExpansions)
223  {
224  returnval->ReadExpansions(infilename);
225  }
226  }
227  return returnval;
228  }
229 
230 
231 
232  /**
233  *
234  */
235  void MeshGraph::ReadGeometry(const std::string& infilename)
236  {
237  TiXmlDocument doc(infilename);
238  bool loadOkay = doc.LoadFile();
239 
240  std::stringstream errstr;
241  errstr << "Unable to load file: " << infilename << " (";
242  errstr << doc.ErrorDesc() << ", line " << doc.ErrorRow()
243  << ", column " << doc.ErrorCol() << ")";
244  ASSERTL0(loadOkay, errstr.str());
245 
246  ReadGeometry(doc);
247  }
248 
249 
250  /**
251  *
252  */
253  void MeshGraph::ReadGeometry(TiXmlDocument &doc)
254  {
255  TiXmlHandle docHandle(&doc);
256  TiXmlElement* mesh = NULL;
257  TiXmlElement* master = NULL; // Master tag within which all data is contained.
258 
259  int err; /// Error value returned by TinyXML.
260 
261  master = doc.FirstChildElement("NEKTAR");
262  ASSERTL0(master, "Unable to find NEKTAR tag in file.");
263 
264  // Find the Mesh tag and same the dim and space attributes
265  mesh = master->FirstChildElement("GEOMETRY");
266 
267  ASSERTL0(mesh, "Unable to find GEOMETRY tag in file.");
268  TiXmlAttribute *attr = mesh->FirstAttribute();
269 
270  // Initialize the mesh and space dimensions to 3 dimensions.
271  // We want to do this each time we read a file, so it should
272  // be done here and not just during class initialization.
273  m_meshPartitioned = false;
274  m_meshDimension = 3;
275  m_spaceDimension = 3;
276 
277  while (attr)
278  {
279  std::string attrName(attr->Name());
280  if (attrName == "DIM")
281  {
282  err = attr->QueryIntValue(&m_meshDimension);
283  ASSERTL1(err==TIXML_SUCCESS, "Unable to read mesh dimension.");
284  }
285  else if (attrName == "SPACE")
286  {
287  err = attr->QueryIntValue(&m_spaceDimension);
288  ASSERTL1(err==TIXML_SUCCESS, "Unable to read space dimension.");
289  }
290  else if (attrName == "PARTITION")
291  {
292  err = attr->QueryIntValue(&m_partition);
293  ASSERTL1(err==TIXML_SUCCESS, "Unable to read partition.");
294  m_meshPartitioned = true;
295  }
296  else
297  {
298  std::string errstr("Unknown attribute: ");
299  errstr += attrName;
300  ASSERTL1(false, errstr.c_str());
301  }
302 
303  // Get the next attribute.
304  attr = attr->Next();
305  }
306 
307  ASSERTL1(m_meshDimension<=m_spaceDimension, "Mesh dimension greater than space dimension");
308 
309  // Now read the vertices
310  TiXmlElement* element = mesh->FirstChildElement("VERTEX");
311  ASSERTL0(element, "Unable to find mesh VERTEX tag in file.");
312 
313  NekDouble xscale,yscale,zscale;
314 
315  // check to see if any scaling parameters are in
316  // attributes and determine these values
317  LibUtilities::AnalyticExpressionEvaluator expEvaluator;
318  //LibUtilities::ExpressionEvaluator expEvaluator;
319  const char *xscal = element->Attribute("XSCALE");
320  if(!xscal)
321  {
322  xscale = 1.0;
323  }
324  else
325  {
326  std::string xscalstr = xscal;
327  int expr_id = expEvaluator.DefineFunction("",xscalstr);
328  xscale = expEvaluator.Evaluate(expr_id);
329  }
330 
331  const char *yscal = element->Attribute("YSCALE");
332  if(!yscal)
333  {
334  yscale = 1.0;
335  }
336  else
337  {
338  std::string yscalstr = yscal;
339  int expr_id = expEvaluator.DefineFunction("",yscalstr);
340  yscale = expEvaluator.Evaluate(expr_id);
341  }
342 
343  const char *zscal = element->Attribute("ZSCALE");
344  if(!zscal)
345  {
346  zscale = 1.0;
347  }
348  else
349  {
350  std::string zscalstr = zscal;
351  int expr_id = expEvaluator.DefineFunction("",zscalstr);
352  zscale = expEvaluator.Evaluate(expr_id);
353  }
354 
355 
356  NekDouble xmove,ymove,zmove;
357 
358  // check to see if any moving parameters are in
359  // attributes and determine these values
360 
361  //LibUtilities::ExpressionEvaluator expEvaluator;
362  const char *xmov = element->Attribute("XMOVE");
363  if(!xmov)
364  {
365  xmove = 0.0;
366  }
367  else
368  {
369  std::string xmovstr = xmov;
370  int expr_id = expEvaluator.DefineFunction("",xmovstr);
371  xmove = expEvaluator.Evaluate(expr_id);
372  }
373 
374  const char *ymov = element->Attribute("YMOVE");
375  if(!ymov)
376  {
377  ymove = 0.0;
378  }
379  else
380  {
381  std::string ymovstr = ymov;
382  int expr_id = expEvaluator.DefineFunction("",ymovstr);
383  ymove = expEvaluator.Evaluate(expr_id);
384  }
385 
386  const char *zmov = element->Attribute("ZMOVE");
387  if(!zmov)
388  {
389  zmove = 0.0;
390  }
391  else
392  {
393  std::string zmovstr = zmov;
394  int expr_id = expEvaluator.DefineFunction("",zmovstr);
395  zmove = expEvaluator.Evaluate(expr_id);
396  }
397 
398  string IsCompressed;
399  element->QueryStringAttribute("COMPRESSED",&IsCompressed);
400 
401  if(IsCompressed.size())
402  {
403  if(boost::iequals(IsCompressed,
404  LibUtilities::CompressData::GetCompressString()))
405  {
406  // Extract the vertex body
407  TiXmlNode* vertexChild = element->FirstChild();
408  ASSERTL0(vertexChild,
409  "Unable to extract the data from the compressed "
410  "vertex tag.");
411 
412  std::string vertexStr;
413  if (vertexChild->Type() == TiXmlNode::TINYXML_TEXT)
414  {
415  vertexStr += vertexChild->ToText()->ValueStr();
416  }
417 
418  std::vector<LibUtilities::MeshVertex> vertData;
419  LibUtilities::CompressData::ZlibDecodeFromBase64Str(
420  vertexStr,vertData);
421 
422  int indx;
423  NekDouble xval, yval, zval;
424  for(int i = 0; i < vertData.size(); ++i)
425  {
426  indx = vertData[i].id;
427  xval = vertData[i].x;
428  yval = vertData[i].y;
429  zval = vertData[i].z;
430 
431  xval = xval*xscale + xmove;
432  yval = yval*yscale + ymove;
433  zval = zval*zscale + zmove;
434 
435  PointGeomSharedPtr vert(
436  MemoryManager<PointGeom>::AllocateSharedPtr(
437  m_spaceDimension, indx, xval, yval, zval));
438 
439  vert->SetGlobalID(indx);
440  m_vertSet[indx] = vert;
441  }
442  }
443  else
444  {
445  ASSERTL0(false,"Compressed formats do not match. Expected :"
446  + LibUtilities::CompressData::GetCompressString()
447  + " but got " + std::string(IsCompressed));
448  }
449  }
450  else
451  {
452  TiXmlElement *vertex = element->FirstChildElement("V");
453 
454  int indx;
455  int nextVertexNumber = -1;
456 
457  while (vertex)
458  {
459  nextVertexNumber++;
460 
461  TiXmlAttribute *vertexAttr = vertex->FirstAttribute();
462  std::string attrName(vertexAttr->Name());
463 
464  ASSERTL0(attrName == "ID", (std::string("Unknown attribute name: ") + attrName).c_str());
465 
466  err = vertexAttr->QueryIntValue(&indx);
467  ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
468 
469  // Now read body of vertex
470  std::string vertexBodyStr;
471 
472  TiXmlNode *vertexBody = vertex->FirstChild();
473 
474  while (vertexBody)
475  {
476  // Accumulate all non-comment body data.
477  if (vertexBody->Type() == TiXmlNode::TINYXML_TEXT)
478  {
479  vertexBodyStr += vertexBody->ToText()->Value();
480  vertexBodyStr += " ";
481  }
482 
483  vertexBody = vertexBody->NextSibling();
484  }
485 
486  ASSERTL0(!vertexBodyStr.empty(), "Vertex definitions must contain vertex data.");
487 
488  // Get vertex data from the data string.
489  NekDouble xval, yval, zval;
490  std::istringstream vertexDataStrm(vertexBodyStr.c_str());
491 
492  try
493  {
494  while(!vertexDataStrm.fail())
495  {
496  vertexDataStrm >> xval >> yval >> zval;
497 
498  xval = xval*xscale + xmove;
499  yval = yval*yscale + ymove;
500  zval = zval*zscale + zmove;
501 
502  // Need to check it here because we may not be
503  // good after the read indicating that there
504  // was nothing to read.
505  if (!vertexDataStrm.fail())
506  {
507  PointGeomSharedPtr vert(MemoryManager<PointGeom>::AllocateSharedPtr(m_spaceDimension, indx, xval, yval, zval));
508  vert->SetGlobalID(indx);
509  m_vertSet[indx] = vert;
510  }
511  }
512  }
513  catch(...)
514  {
515  ASSERTL0(false, "Unable to read VERTEX data.");
516  }
517 
518  vertex = vertex->NextSiblingElement("V");
519  }
520  }
521  }
522 
523 
524  /**
525  * Read the geometry-related information from the given file. This
526  * information is located within the XML tree under
527  * <NEKTAR><GEOMETRY><GEOMINFO>.
528  * @param infilename Filename of XML file.
529  */
530  void MeshGraph::ReadGeometryInfo(const std::string &infilename)
531  {
532  TiXmlDocument doc(infilename);
533  bool loadOkay = doc.LoadFile();
534 
535  std::stringstream errstr;
536  errstr << "Unable to load file: " << infilename << std::endl;
537  errstr << "Reason: " << doc.ErrorDesc() << std::endl;
538  errstr << "Position: Line " << doc.ErrorRow() << ", Column " << doc.ErrorCol() << std::endl;
539  ASSERTL0(loadOkay, errstr.str());
540 
541  ReadGeometryInfo(doc);
542  }
543 
544 
545  /**
546  * Read the geometry-related information from the given XML document.
547  * This information is located within the XML tree under
548  * <NEKTAR><GEOMETRY><GEOMINFO>.
549  * @param doc XML document.
550  */
551  void MeshGraph::ReadGeometryInfo(TiXmlDocument &doc)
552  {
553  TiXmlElement *master = doc.FirstChildElement("NEKTAR");
554  ASSERTL0(master, "Unable to find NEKTAR tag in file.");
555 
556  // Find the Expansions tag
557  TiXmlElement *geomTag = master->FirstChildElement("GEOMETRY");
558  ASSERTL0(geomTag, "Unable to find GEOMETRY tag in file.");
559 
560  // See if we have GEOMINFO. If there is none, it's fine.
561  TiXmlElement *geomInfoTag = geomTag->FirstChildElement("GEOMINFO");
562  if (!geomInfoTag) return;
563 
564  TiXmlElement *infoItem = geomInfoTag->FirstChildElement("I");
565 
566  // Multiple nodes will only occur if there is a comment in between
567  // definitions.
568  while (infoItem)
569  {
570  std::string geomProperty = infoItem->Attribute("PROPERTY");
571  std::string geomValue = infoItem->Attribute("VALUE");
572  GeomInfoMap::iterator x = m_geomInfo.find(geomProperty);
573 
574  ASSERTL0(x == m_geomInfo.end(),
575  "Property " + geomProperty + " already specified.");
576  m_geomInfo[geomProperty] = geomValue;
577  infoItem = infoItem->NextSiblingElement("I");
578  }
579  }
580 
581 
582  /**
583  *
584  */
585  void MeshGraph::ReadExpansions(const std::string& infilename)
586  {
587  TiXmlDocument doc(infilename);
588  bool loadOkay = doc.LoadFile();
589 
590  std::stringstream errstr;
591  errstr << "Unable to load file: " << infilename << std::endl;
592  errstr << "Reason: " << doc.ErrorDesc() << std::endl;
593  errstr << "Position: Line " << doc.ErrorRow() << ", Column " << doc.ErrorCol() << std::endl;
594  ASSERTL0(loadOkay, errstr.str());
595 
596  ReadExpansions(doc);
597  }
598 
599 
600  /**
601  *
602  */
603  void MeshGraph::ReadExpansions(TiXmlDocument &doc)
604  {
605  TiXmlElement *master = doc.FirstChildElement("NEKTAR");
606  ASSERTL0(master, "Unable to find NEKTAR tag in file.");
607 
608  // Find the Expansions tag
609  TiXmlElement *expansionTypes = master->FirstChildElement("EXPANSIONS");
610  ASSERTL0(expansionTypes, "Unable to find EXPANSIONS tag in file.");
611 
612  if(expansionTypes)
613  {
614  // Find the Expansion type
615  TiXmlElement *expansion = expansionTypes->FirstChildElement();
616  std::string expType = expansion->Value();
617 
618  if(expType == "E")
619  {
620  int i;
621  ExpansionMapShPtr expansionMap;
622 
623  /// Expansiontypes will contain composite,
624  /// nummodes, and expansiontype (eModified, or
625  /// eOrthogonal) Or a full list of data of
626  /// basistype, nummodes, pointstype, numpoints;
627 
628  /// Expansiontypes may also contain a list of
629  /// fields that this expansion relates to. If this
630  /// does not exist the variable is only set to
631  /// "DefaultVar".
632 
633  while (expansion)
634  {
635 
636  const char *fStr = expansion->Attribute("FIELDS");
637  std::vector<std::string> fieldStrings;
638 
639  if(fStr) // extract other fields.
640  {
641  std::string fieldStr = fStr;
642  bool valid = ParseUtils::GenerateOrderedStringVector(fieldStr.c_str(),fieldStrings);
643  ASSERTL0(valid,"Unable to correctly parse the field string in ExpansionTypes.");
644  }
645 
646  // check to see if m_expasionVectorShPtrMap has
647  // already been intiailised and if not intiailse
648  // vector.
649  if(m_expansionMapShPtrMap.count("DefaultVar") == 0) // no previous definitions
650  {
651  expansionMap = SetUpExpansionMap();
652 
653  m_expansionMapShPtrMap["DefaultVar"] = expansionMap;
654 
655  // make sure all fields in this search point
656  // to same expansion vector;
657  for(i = 0; i < fieldStrings.size(); ++i)
658  {
659  m_expansionMapShPtrMap[fieldStrings[i]] = expansionMap;
660  }
661  }
662  else // default variable is defined
663  {
664 
665  if(fieldStrings.size()) // fields are defined
666  {
667  //see if field exists
668  if(m_expansionMapShPtrMap.count(fieldStrings[0]))
669  {
670  expansionMap = m_expansionMapShPtrMap.find(fieldStrings[0])->second;
671  }
672  else
673  {
674  expansionMap = SetUpExpansionMap();
675  // make sure all fields in this search point
676  // to same expansion vector;
677  for(i = 0; i < fieldStrings.size(); ++i)
678  {
679  if(m_expansionMapShPtrMap.count(fieldStrings[i]) == 0)
680  {
681  m_expansionMapShPtrMap[fieldStrings[i]] = expansionMap;
682  }
683  else
684  {
685  ASSERTL0(false,"Expansion vector for this field is already setup");
686  }
687  }
688  }
689  }
690  else // use default variable list
691  {
692  expansionMap = m_expansionMapShPtrMap.find("DefaultVar")->second;
693  }
694 
695  }
696 
697  /// Mandatory components...optional are to follow later.
698  std::string compositeStr = expansion->Attribute("COMPOSITE");
699  ASSERTL0(compositeStr.length() > 3, "COMPOSITE must be specified in expansion definition");
700  int beg = compositeStr.find_first_of("[");
701  int end = compositeStr.find_first_of("]");
702  std::string compositeListStr = compositeStr.substr(beg+1,end-beg-1);
703 
704  CompositeMap compositeVector;
705  GetCompositeList(compositeListStr, compositeVector);
706 
707  bool useExpansionType = false;
708  ExpansionType expansion_type;
709  int num_modes;
710 
711  LibUtilities::BasisKeyVector basiskeyvec;
712  const char * tStr = expansion->Attribute("TYPE");
713 
714  if(tStr) // use type string to define expansion
715  {
716  std::string typeStr = tStr;
717  const std::string* begStr = kExpansionTypeStr;
718  const std::string* endStr = kExpansionTypeStr+eExpansionTypeSize;
719  const std::string* expStr = std::find(begStr, endStr, typeStr);
720 
721  ASSERTL0(expStr != endStr, "Invalid expansion type.");
722  expansion_type = (ExpansionType)(expStr - begStr);
723 
724 
725  /// \todo solvers break the pattern 'instantiate Session -> instantiate MeshGraph'
726  /// and parse command line arguments by themselves; one needs to unify command
727  /// line arguments handling.
728  /// Solvers tend to call MeshGraph::Read statically -> m_session
729  /// is not defined -> no info about command line arguments presented
730  /// ASSERTL0(m_session != 0, "One needs to instantiate SessionReader first");
731 
732  const char *nStr = expansion->Attribute("NUMMODES");
733  ASSERTL0(nStr,"NUMMODES was not defined in EXPANSION section of input");
734  std::string nummodesStr = nStr;
735 
736  // ASSERTL0(m_session,"Session should be defined to evaluate nummodes ");
737  if (m_session)
738  {
739  LibUtilities::Equation nummodesEqn(m_session, nummodesStr);
740  num_modes = (int) nummodesEqn.Evaluate();
741  }
742  else
743  {
744  num_modes = boost::lexical_cast<int>(nummodesStr);
745  }
746 
747  useExpansionType = true;
748  }
749  else // assume expansion is defined individually
750  {
751  // Extract the attributes.
752  const char *bTypeStr = expansion->Attribute("BASISTYPE");
753  ASSERTL0(bTypeStr,"TYPE or BASISTYPE was not defined in EXPANSION section of input");
754  std::string basisTypeStr = bTypeStr;
755 
756  // interpret the basis type string.
757  std::vector<std::string> basisStrings;
758  std::vector<LibUtilities::BasisType> basis;
759  bool valid = ParseUtils::GenerateOrderedStringVector(basisTypeStr.c_str(), basisStrings);
760  ASSERTL0(valid, "Unable to correctly parse the basis types.");
761  for (vector<std::string>::size_type i = 0; i < basisStrings.size(); i++)
762  {
763  valid = false;
764  for (unsigned int j = 0; j < LibUtilities::SIZE_BasisType; j++)
765  {
766  if (LibUtilities::BasisTypeMap[j] == basisStrings[i])
767  {
768  basis.push_back((LibUtilities::BasisType) j);
769  valid = true;
770  break;
771  }
772  }
773  ASSERTL0(valid, std::string("Unable to correctly parse the basis type: ").append(basisStrings[i]).c_str());
774  }
775  const char *nModesStr = expansion->Attribute("NUMMODES");
776  ASSERTL0(nModesStr,"NUMMODES was not defined in EXPANSION section of input");
777 
778  std::string numModesStr = nModesStr;
779  std::vector<unsigned int> numModes;
780  valid = ParseUtils::GenerateOrderedVector(numModesStr.c_str(), numModes);
781  ASSERTL0(valid, "Unable to correctly parse the number of modes.");
782  ASSERTL0(numModes.size() == basis.size(),"information for num modes does not match the number of basis");
783 
784  const char *pTypeStr = expansion->Attribute("POINTSTYPE");
785  ASSERTL0(pTypeStr,"POINTSTYPE was not defined in EXPANSION section of input");
786  std::string pointsTypeStr = pTypeStr;
787  // interpret the points type string.
788  std::vector<std::string> pointsStrings;
789  std::vector<LibUtilities::PointsType> points;
790  valid = ParseUtils::GenerateOrderedStringVector(pointsTypeStr.c_str(), pointsStrings);
791  ASSERTL0(valid, "Unable to correctly parse the points types.");
792  for (vector<std::string>::size_type i = 0; i < pointsStrings.size(); i++)
793  {
794  valid = false;
795  for (unsigned int j = 0; j < LibUtilities::SIZE_PointsType; j++)
796  {
797  if (LibUtilities::kPointsTypeStr[j] == pointsStrings[i])
798  {
799  points.push_back((LibUtilities::PointsType) j);
800  valid = true;
801  break;
802  }
803  }
804  ASSERTL0(valid, std::string("Unable to correctly parse the points type: ").append(pointsStrings[i]).c_str());
805  }
806 
807  const char *nPointsStr = expansion->Attribute("NUMPOINTS");
808  ASSERTL0(nPointsStr,"NUMPOINTS was not defined in EXPANSION section of input");
809  std::string numPointsStr = nPointsStr;
810  std::vector<unsigned int> numPoints;
811  valid = ParseUtils::GenerateOrderedVector(numPointsStr.c_str(), numPoints);
812  ASSERTL0(valid, "Unable to correctly parse the number of points.");
813  ASSERTL0(numPoints.size() == numPoints.size(),"information for num points does not match the number of basis");
814 
815  for(int i = 0; i < basis.size(); ++i)
816  {
817  //Generate Basis key using information
818  const LibUtilities::PointsKey pkey(numPoints[i],points[i]);
819  basiskeyvec.push_back(LibUtilities::BasisKey(basis[i],numModes[i],pkey));
820  }
821  }
822 
823  // Now have composite and basiskeys. Cycle through
824  // all composites for the geomShPtrs and set the modes
825  // and types for the elements contained in the element
826  // list.
827  CompositeMapIter compVecIter;
828  for (compVecIter = compositeVector.begin(); compVecIter != compositeVector.end(); ++compVecIter)
829  {
830  GeometryVectorIter geomVecIter;
831  for (geomVecIter = (compVecIter->second)->begin(); geomVecIter != (compVecIter->second)->end(); ++geomVecIter)
832  {
833  ExpansionMapIter x = expansionMap->find((*geomVecIter)->GetGlobalID());
834  ASSERTL0(x != expansionMap->end(), "Expansion not found!!");
835  if(useExpansionType)
836  {
837  (x->second)->m_basisKeyVector = MeshGraph::DefineBasisKeyFromExpansionType(*geomVecIter,expansion_type,num_modes);
838  }
839  else
840  {
841  ASSERTL0((*geomVecIter)->GetShapeDim() == basiskeyvec.size()," There is an incompatible expansion dimension with geometry dimension");
842  (x->second)->m_basisKeyVector = basiskeyvec;
843  }
844  }
845  }
846 
847  expansion = expansion->NextSiblingElement("E");
848  }
849  }
850  else if(expType == "H")
851  {
852  int i;
853  ExpansionMapShPtr expansionMap;
854 
855  while (expansion)
856  {
857 
858  const char *fStr = expansion->Attribute("FIELDS");
859  std::vector<std::string> fieldStrings;
860 
861  if(fStr) // extract other fields.
862  {
863  std::string fieldStr = fStr;
864  bool valid = ParseUtils::GenerateOrderedStringVector(fieldStr.c_str(),fieldStrings);
865  ASSERTL0(valid,"Unable to correctly parse the field string in ExpansionTypes.");
866  }
867 
868  // check to see if m_expasionVectorShPtrMap has
869  // already been intiailised and if not intiailse
870  // vector.
871  if(m_expansionMapShPtrMap.count("DefaultVar") == 0) // no previous definitions
872  {
873  expansionMap = SetUpExpansionMap();
874 
875  m_expansionMapShPtrMap["DefaultVar"] = expansionMap;
876 
877  // make sure all fields in this search point
878  // to same expansion vector;
879  for(i = 0; i < fieldStrings.size(); ++i)
880  {
881  m_expansionMapShPtrMap[fieldStrings[i]] = expansionMap;
882  }
883  }
884  else // default variable is defined
885  {
886 
887  if(fieldStrings.size()) // fields are defined
888  {
889  //see if field exists
890  if(m_expansionMapShPtrMap.count(fieldStrings[0]))
891  {
892  expansionMap = m_expansionMapShPtrMap.find(fieldStrings[0])->second;
893  }
894  else
895  {
896  expansionMap = SetUpExpansionMap();
897  // make sure all fields in this search point
898  // to same expansion vector;
899  for(i = 0; i < fieldStrings.size(); ++i)
900  {
901  if(m_expansionMapShPtrMap.count(fieldStrings[i]) == 0)
902  {
903  m_expansionMapShPtrMap[fieldStrings[i]] = expansionMap;
904  }
905  else
906  {
907  ASSERTL0(false,"Expansion vector for this field is already setup");
908  }
909  }
910  }
911  }
912  else // use default variable list
913  {
914  expansionMap = m_expansionMapShPtrMap.find("DefaultVar")->second;
915  }
916 
917  }
918 
919  /// Mandatory components...optional are to follow later.
920  std::string compositeStr = expansion->Attribute("COMPOSITE");
921  ASSERTL0(compositeStr.length() > 3, "COMPOSITE must be specified in expansion definition");
922  int beg = compositeStr.find_first_of("[");
923  int end = compositeStr.find_first_of("]");
924  std::string compositeListStr = compositeStr.substr(beg+1,end-beg-1);
925 
926  CompositeMap compositeVector;
927  GetCompositeList(compositeListStr, compositeVector);
928 
929  ExpansionType expansion_type_x = eNoExpansionType;
930  ExpansionType expansion_type_y = eNoExpansionType;
931  ExpansionType expansion_type_z = eNoExpansionType;
932  int num_modes_x = 0;
933  int num_modes_y = 0;
934  int num_modes_z = 0;
935 
936  LibUtilities::BasisKeyVector basiskeyvec;
937 
938  const char * tStr_x = expansion->Attribute("TYPE-X");
939 
940  if(tStr_x) // use type string to define expansion
941  {
942  std::string typeStr = tStr_x;
943  const std::string* begStr = kExpansionTypeStr;
944  const std::string* endStr = kExpansionTypeStr+eExpansionTypeSize;
945  const std::string* expStr = std::find(begStr, endStr, typeStr);
946 
947  ASSERTL0(expStr != endStr, "Invalid expansion type.");
948  expansion_type_x = (ExpansionType)(expStr - begStr);
949 
950  const char *nStr = expansion->Attribute("NUMMODES-X");
951  ASSERTL0(nStr,"NUMMODES-X was not defined in EXPANSION section of input");
952  std::string nummodesStr = nStr;
953 
954  // ASSERTL0(m_session,"Session should be defined to evaluate nummodes ");
955 
956  if (m_session)
957  {
958  LibUtilities::Equation nummodesEqn(m_session, nummodesStr);
959  num_modes_x = (int) nummodesEqn.Evaluate();
960  }
961  else
962  {
963  num_modes_x = boost::lexical_cast<int>(nummodesStr);
964  }
965 
966  }
967 
968  const char * tStr_y = expansion->Attribute("TYPE-Y");
969 
970  if(tStr_y) // use type string to define expansion
971  {
972  std::string typeStr = tStr_y;
973  const std::string* begStr = kExpansionTypeStr;
974  const std::string* endStr = kExpansionTypeStr+eExpansionTypeSize;
975  const std::string* expStr = std::find(begStr, endStr, typeStr);
976 
977  ASSERTL0(expStr != endStr, "Invalid expansion type.");
978  expansion_type_y = (ExpansionType)(expStr - begStr);
979 
980  const char *nStr = expansion->Attribute("NUMMODES-Y");
981  ASSERTL0(nStr,"NUMMODES-Y was not defined in EXPANSION section of input");
982  std::string nummodesStr = nStr;
983 
984  // ASSERTL0(m_session,"Session should be defined to evaluate nummodes ");
985  if (m_session)
986  {
987  LibUtilities::Equation nummodesEqn(m_session, nummodesStr);
988  num_modes_y = (int) nummodesEqn.Evaluate();
989  }
990  else
991  {
992  num_modes_y = boost::lexical_cast<int>(nummodesStr);
993  }
994 
995  }
996 
997  const char * tStr_z = expansion->Attribute("TYPE-Z");
998 
999  if(tStr_z) // use type string to define expansion
1000  {
1001  std::string typeStr = tStr_z;
1002  const std::string* begStr = kExpansionTypeStr;
1003  const std::string* endStr = kExpansionTypeStr+eExpansionTypeSize;
1004  const std::string* expStr = std::find(begStr, endStr, typeStr);
1005 
1006  ASSERTL0(expStr != endStr, "Invalid expansion type.");
1007  expansion_type_z = (ExpansionType)(expStr - begStr);
1008 
1009  const char *nStr = expansion->Attribute("NUMMODES-Z");
1010  ASSERTL0(nStr,"NUMMODES-Z was not defined in EXPANSION section of input");
1011  std::string nummodesStr = nStr;
1012 
1013  // ASSERTL0(m_session,"Session should be defined to evaluate nummodes ");
1014  if (m_session)
1015  {
1016  LibUtilities::Equation nummodesEqn(m_session, nummodesStr);
1017  num_modes_z = (int) nummodesEqn.Evaluate();
1018  }
1019  else
1020  {
1021  num_modes_z = boost::lexical_cast<int>(nummodesStr);
1022  }
1023 
1024  }
1025 
1026  CompositeMapIter compVecIter;
1027  for (compVecIter = compositeVector.begin(); compVecIter != compositeVector.end(); ++compVecIter)
1028  {
1029  GeometryVectorIter geomVecIter;
1030  for (geomVecIter = (compVecIter->second)->begin(); geomVecIter != (compVecIter->second)->end(); ++geomVecIter)
1031  {
1032  ExpansionMapIter expVecIter;
1033  for (expVecIter = expansionMap->begin(); expVecIter != expansionMap->end(); ++expVecIter)
1034  {
1035 
1036  (expVecIter->second)->m_basisKeyVector = DefineBasisKeyFromExpansionTypeHomo(*geomVecIter,
1037  expansion_type_x,
1038  expansion_type_y,
1039  expansion_type_z,
1040  num_modes_x,
1041  num_modes_y,
1042  num_modes_z);
1043  }
1044  }
1045  }
1046 
1047  expansion = expansion->NextSiblingElement("H");
1048  }
1049  }
1050  else if(expType == "ELEMENTS") // Reading a file with the expansion definition
1051  {
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);
1057  cout << " Number of elements: " << fielddefs.size() << endl;
1058  SetExpansions(fielddefs);
1059  }
1060  else if(expType == "F")
1061  {
1062  ASSERTL0(expansion->Attribute("FILE"),
1063  "Attribute FILE expected for type F expansion");
1064  std::string filenameStr = expansion->Attribute("FILE");
1065  ASSERTL0(!filenameStr.empty(),
1066  "A filename must be specified for the FILE "
1067  "attribute of expansion");
1068 
1069  std::vector<LibUtilities::FieldDefinitionsSharedPtr> fielddefs;
1070  LibUtilities::FieldIOSharedPtr f =
1071  LibUtilities::FieldIO::CreateForFile(
1072  m_session, filenameStr);
1073  f->Import(filenameStr, fielddefs);
1074  SetExpansions(fielddefs);
1075  }
1076  else
1077  {
1078  ASSERTL0(false,"Expansion type not defined");
1079  }
1080  }
1081  }
1082 
1083 
1084  /**
1085  *
1086  */
1087  void MeshGraph::ReadDomain(TiXmlDocument &doc)
1088  {
1089  TiXmlHandle docHandle(&doc);
1090 
1091  TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
1092  TiXmlElement* domain = NULL;
1093 
1094  ASSERTL0(mesh, "Unable to find GEOMETRY tag in file.");
1095 
1096  /// Look for data in DOMAIN block.
1097  domain = mesh->FirstChildElement("DOMAIN");
1098 
1099  ASSERTL0(domain, "Unable to find DOMAIN tag in file.");
1100 
1101  /// Elements are of the form: "<D ID = "N"> ... </D>".
1102  /// Read the ID field first.
1103  TiXmlElement *multidomains = domain->FirstChildElement("D");
1104 
1105  if(multidomains)
1106  {
1107  int nextDomainNumber = 0;
1108  while (multidomains)
1109  {
1110  int indx;
1111  int err = multidomains->QueryIntAttribute("ID", &indx);
1112  ASSERTL0(err == TIXML_SUCCESS,
1113  "Unable to read attribute ID in Domain.");
1114 
1115 
1116  TiXmlNode* elementChild = multidomains->FirstChild();
1117  while(elementChild && elementChild->Type() != TiXmlNode::TINYXML_TEXT)
1118  {
1119  elementChild = elementChild->NextSibling();
1120  }
1121 
1122  ASSERTL0(elementChild, "Unable to read DOMAIN body.");
1123  std::string elementStr = elementChild->ToText()->ValueStr();
1124 
1125  elementStr = elementStr.substr(elementStr.find_first_not_of(" "));
1126 
1127  std::string::size_type indxBeg = elementStr.find_first_of('[') + 1;
1128  std::string::size_type indxEnd = elementStr.find_last_of(']') - 1;
1129  std::string indxStr = elementStr.substr(indxBeg, indxEnd - indxBeg + 1);
1130 
1131  ASSERTL0(!indxStr.empty(), "Unable to read domain's composite index (index missing?).");
1132 
1133  // Read the domain composites.
1134  // Parse the composites into a list.
1135  CompositeMap unrollDomain;
1136  GetCompositeList(indxStr, unrollDomain);
1137  m_domain.push_back(unrollDomain);
1138 
1139  ASSERTL0(!m_domain[nextDomainNumber++].empty(), (std::string("Unable to obtain domain's referenced composite: ") + indxStr).c_str());
1140 
1141  /// Keep looking
1142  multidomains = multidomains->NextSiblingElement("D");
1143  }
1144 
1145  }
1146  else // previous definition of just one composite
1147  {
1148 
1149  // find the non comment portion of the body.
1150  TiXmlNode* elementChild = domain->FirstChild();
1151  while(elementChild && elementChild->Type() != TiXmlNode::TINYXML_TEXT)
1152  {
1153  elementChild = elementChild->NextSibling();
1154  }
1155 
1156  ASSERTL0(elementChild, "Unable to read DOMAIN body.");
1157  std::string elementStr = elementChild->ToText()->ValueStr();
1158 
1159  elementStr = elementStr.substr(elementStr.find_first_not_of(" "));
1160 
1161  std::string::size_type indxBeg = elementStr.find_first_of('[') + 1;
1162  std::string::size_type indxEnd = elementStr.find_last_of(']') - 1;
1163  std::string indxStr = elementStr.substr(indxBeg, indxEnd - indxBeg + 1);
1164 
1165  ASSERTL0(!indxStr.empty(), "Unable to read domain's composite index (index missing?).");
1166 
1167  // Read the domain composites.
1168  // Parse the composites into a list.
1169  CompositeMap fullDomain;
1170  GetCompositeList(indxStr, fullDomain);
1171  m_domain.push_back(fullDomain);
1172 
1173  ASSERTL0(!m_domain[0].empty(), (std::string("Unable to obtain domain's referenced composite: ") + indxStr).c_str());
1174  }
1175  }
1176 
1177 
1178  /**
1179  *
1180  */
1181  void MeshGraph::ReadCurves(TiXmlDocument &doc)
1182  {
1183  /// We know we have it since we made it this far.
1184  TiXmlHandle docHandle(&doc);
1185  TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
1186  TiXmlElement* field = NULL;
1187 
1188  // check to see if any scaling parameters are in
1189  // attributes and determine these values
1190  TiXmlElement* element = mesh->FirstChildElement("VERTEX");
1191  ASSERTL0(element, "Unable to find mesh VERTEX tag in file.");
1192 
1193  NekDouble xscale,yscale,zscale;
1194 
1195  LibUtilities::AnalyticExpressionEvaluator expEvaluator;
1196  const char *xscal = element->Attribute("XSCALE");
1197  if(!xscal)
1198  {
1199  xscale = 1.0;
1200  }
1201  else
1202  {
1203  std::string xscalstr = xscal;
1204  int expr_id = expEvaluator.DefineFunction("",xscalstr);
1205  xscale = expEvaluator.Evaluate(expr_id);
1206  }
1207 
1208  const char *yscal = element->Attribute("YSCALE");
1209  if(!yscal)
1210  {
1211  yscale = 1.0;
1212  }
1213  else
1214  {
1215  std::string yscalstr = yscal;
1216  int expr_id = expEvaluator.DefineFunction("",yscalstr);
1217  yscale = expEvaluator.Evaluate(expr_id);
1218  }
1219 
1220  const char *zscal = element->Attribute("ZSCALE");
1221  if(!zscal)
1222  {
1223  zscale = 1.0;
1224  }
1225  else
1226  {
1227  std::string zscalstr = zscal;
1228  int expr_id = expEvaluator.DefineFunction("",zscalstr);
1229  zscale = expEvaluator.Evaluate(expr_id);
1230  }
1231 
1232  NekDouble xmove,ymove,zmove;
1233 
1234  // check to see if any moving parameters are in
1235  // attributes and determine these values
1236 
1237  //LibUtilities::ExpressionEvaluator expEvaluator;
1238  const char *xmov = element->Attribute("XMOVE");
1239  if(!xmov)
1240  {
1241  xmove = 0.0;
1242  }
1243  else
1244  {
1245  std::string xmovstr = xmov;
1246  int expr_id = expEvaluator.DefineFunction("",xmovstr);
1247  xmove = expEvaluator.Evaluate(expr_id);
1248  }
1249 
1250  const char *ymov = element->Attribute("YMOVE");
1251  if(!ymov)
1252  {
1253  ymove = 0.0;
1254  }
1255  else
1256  {
1257  std::string ymovstr = ymov;
1258  int expr_id = expEvaluator.DefineFunction("",ymovstr);
1259  ymove = expEvaluator.Evaluate(expr_id);
1260  }
1261 
1262  const char *zmov = element->Attribute("ZMOVE");
1263  if(!zmov)
1264  {
1265  zmove = 0.0;
1266  }
1267  else
1268  {
1269  std::string zmovstr = zmov;
1270  int expr_id = expEvaluator.DefineFunction("",zmovstr);
1271  zmove = expEvaluator.Evaluate(expr_id);
1272  }
1273 
1274  int err;
1275 
1276  /// Look for elements in CURVE block.
1277  field = mesh->FirstChildElement("CURVED");
1278 
1279  if(!field) //return if no curved entities
1280  {
1281  return;
1282  }
1283 
1284  string IsCompressed;
1285  field->QueryStringAttribute("COMPRESSED",&IsCompressed);
1286 
1287  if(IsCompressed.size())
1288  {
1289  ASSERTL0(boost::iequals(IsCompressed,
1290  LibUtilities::CompressData::GetCompressString()),
1291  "Compressed formats do not match. Expected :"
1292  + LibUtilities::CompressData::GetCompressString()
1293  + " but got "
1294  + boost::lexical_cast<std::string>(IsCompressed));
1295 
1296  std::vector<LibUtilities::MeshCurvedInfo> edginfo;
1297  std::vector<LibUtilities::MeshCurvedInfo> facinfo;
1298  LibUtilities::MeshCurvedPts cpts;
1299 
1300  // read edge, face info and curved poitns.
1301  TiXmlElement *x = field->FirstChildElement();
1302  while(x)
1303  {
1304  const char *entitytype = x->Value();
1305  // read in edge or face info
1306  if(boost::iequals(entitytype,"E"))
1307  {
1308  // read in data
1309  std::string elmtStr;
1310  TiXmlNode* child = x->FirstChild();
1311 
1312  if (child->Type() == TiXmlNode::TINYXML_TEXT)
1313  {
1314  elmtStr += child->ToText()->ValueStr();
1315  }
1316 
1317  LibUtilities::CompressData::ZlibDecodeFromBase64Str(
1318  elmtStr,edginfo);
1319  }
1320  else if(boost::iequals(entitytype,"F"))
1321  {
1322  // read in data
1323  std::string elmtStr;
1324  TiXmlNode* child = x->FirstChild();
1325 
1326  if (child->Type() == TiXmlNode::TINYXML_TEXT)
1327  {
1328  elmtStr += child->ToText()->ValueStr();
1329  }
1330 
1331  LibUtilities::CompressData::ZlibDecodeFromBase64Str(
1332  elmtStr,facinfo);
1333  }
1334  else if(boost::iequals(entitytype,"DATAPOINTS"))
1335  {
1336  NekInt id;
1337  ASSERTL0(x->Attribute("ID", &id),
1338  "Failed to get ID from PTS section");
1339  cpts.id = id;
1340 
1341  // read in data
1342  std::string elmtStr;
1343 
1344  TiXmlElement* DataIdx =
1345  x->FirstChildElement("INDEX");
1346  ASSERTL0(DataIdx,
1347  "Cannot read data index tag in compressed "
1348  "curved section");
1349 
1350  TiXmlNode* child = DataIdx->FirstChild();
1351  if (child->Type() == TiXmlNode::TINYXML_TEXT)
1352  {
1353  elmtStr = child->ToText()->ValueStr();
1354  }
1355 
1356  LibUtilities::CompressData::ZlibDecodeFromBase64Str(
1357  elmtStr,cpts.index);
1358 
1359  TiXmlElement* DataPts =
1360  x->FirstChildElement("POINTS");
1361  ASSERTL0(DataPts,
1362  "Cannot read data pts tag in compressed "
1363  "curved section");
1364 
1365  child = DataPts->FirstChild();
1366  if (child->Type() == TiXmlNode::TINYXML_TEXT)
1367  {
1368  elmtStr = child->ToText()->ValueStr();
1369  }
1370 
1371  LibUtilities::CompressData::ZlibDecodeFromBase64Str(
1372  elmtStr,cpts.pts);
1373  }
1374  else
1375  {
1376  ASSERTL0(false,"Unknown tag in curved section");
1377  }
1378  x = x->NextSiblingElement();
1379  }
1380 
1381  // rescale (x,y,z) points;
1382  for(int i = 0; i > cpts.pts.size(); ++i)
1383  {
1384  cpts.pts[i].x = xscale*cpts.pts[i].x + xmove;
1385  cpts.pts[i].y = yscale*cpts.pts[i].y + ymove;
1386  cpts.pts[i].z = zscale*cpts.pts[i].z + zmove;
1387  }
1388 
1389  for(int i = 0; i < edginfo.size(); ++i)
1390  {
1391  int edgeid = edginfo[i].entityid;
1392  LibUtilities::PointsType ptype;
1393 
1394  CurveSharedPtr curve(
1395  MemoryManager<Curve>::AllocateSharedPtr(
1396  edgeid, ptype = (LibUtilities::PointsType)
1397  edginfo[i].ptype));
1398 
1399  // load points
1400  int offset = edginfo[i].ptoffset;
1401  for(int j = 0; j < edginfo[i].npoints; ++j)
1402  {
1403  int idx = cpts.index[offset+j];
1404 
1405  PointGeomSharedPtr vert(
1406  MemoryManager<PointGeom>::AllocateSharedPtr(
1407  m_meshDimension, edginfo[i].id,
1408  cpts.pts[idx].x, cpts.pts[idx].y,
1409  cpts.pts[idx].z));
1410  curve->m_points.push_back(vert);
1411  }
1412 
1413  m_curvedEdges[edgeid] = curve;
1414  }
1415 
1416  for(int i = 0; i < facinfo.size(); ++i)
1417  {
1418  int faceid = facinfo[i].entityid;
1419  LibUtilities::PointsType ptype;
1420 
1421  CurveSharedPtr curve(
1422  MemoryManager<Curve>::AllocateSharedPtr(
1423  faceid, ptype = (LibUtilities::PointsType)
1424  facinfo[i].ptype));
1425 
1426  int offset = facinfo[i].ptoffset;
1427  for(int j = 0; j < facinfo[i].npoints; ++j)
1428  {
1429  int idx = cpts.index[offset+j];
1430 
1431  PointGeomSharedPtr vert(MemoryManager<PointGeom>::
1432  AllocateSharedPtr(m_meshDimension,
1433  facinfo[i].id,
1434  cpts.pts[idx].x,
1435  cpts.pts[idx].y,
1436  cpts.pts[idx].z));
1437  curve->m_points.push_back(vert);
1438  }
1439 
1440  m_curvedFaces[faceid] = curve;
1441  }
1442  }
1443  else
1444  {
1445  /// All curves are of the form: "<? ID="#" TYPE="GLL OR other
1446  /// points type" NUMPOINTS="#"> ... </?>", with ? being an
1447  /// element type (either E or F).
1448 
1449  TiXmlElement *edgelement = field->FirstChildElement("E");
1450 
1451  int edgeindx, edgeid;
1452  int nextEdgeNumber = -1;
1453 
1454  while(edgelement)
1455  {
1456  /// These should be ordered.
1457  nextEdgeNumber++;
1458 
1459  std::string edge(edgelement->ValueStr());
1460  ASSERTL0(edge == "E", (std::string("Unknown 3D curve type:") + edge).c_str());
1461 
1462  /// Read id attribute.
1463  err = edgelement->QueryIntAttribute("ID", &edgeindx);
1464  ASSERTL0(err == TIXML_SUCCESS, "Unable to read curve attribute ID.");
1465 
1466  /// Read edge id attribute.
1467  err = edgelement->QueryIntAttribute("EDGEID", &edgeid);
1468  ASSERTL0(err == TIXML_SUCCESS, "Unable to read curve attribute EDGEID.");
1469 
1470  /// Read text edgelement description.
1471  std::string elementStr;
1472  TiXmlNode* elementChild = edgelement->FirstChild();
1473 
1474  while(elementChild)
1475  {
1476  // Accumulate all non-comment element data
1477  if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
1478  {
1479  elementStr += elementChild->ToText()->ValueStr();
1480  elementStr += " ";
1481  }
1482  elementChild = elementChild->NextSibling();
1483  }
1484 
1485  ASSERTL0(!elementStr.empty(), "Unable to read curve description body.");
1486 
1487  /// Parse out the element components corresponding to type of element.
1488  if (edge == "E")
1489  {
1490  int numPts=0;
1491  // Determine the points type
1492  std::string typeStr = edgelement->Attribute("TYPE");
1493  ASSERTL0(!typeStr.empty(), "TYPE must be specified in " "points definition");
1494 
1495  LibUtilities::PointsType type;
1496  const std::string* begStr = LibUtilities::kPointsTypeStr;
1497  const std::string* endStr = LibUtilities::kPointsTypeStr + LibUtilities::SIZE_PointsType;
1498  const std::string* ptsStr = std::find(begStr, endStr, typeStr);
1499 
1500  ASSERTL0(ptsStr != endStr, "Invalid points type.");
1501  type = (LibUtilities::PointsType)(ptsStr - begStr);
1502 
1503  //Determine the number of points
1504  err = edgelement->QueryIntAttribute("NUMPOINTS", &numPts);
1505  ASSERTL0(err == TIXML_SUCCESS, "Unable to read curve attribute NUMPOINTS.");
1506  CurveSharedPtr curve(MemoryManager<Curve>::AllocateSharedPtr(edgeid, type));
1507 
1508  // Read points (x, y, z)
1509  NekDouble xval, yval, zval;
1510  std::istringstream elementDataStrm(elementStr.c_str());
1511  try
1512  {
1513  while(!elementDataStrm.fail())
1514  {
1515  elementDataStrm >> xval >> yval >> zval;
1516 
1517  xval = xval*xscale + xmove;
1518  yval = yval*yscale + ymove;
1519  zval = zval*zscale + zmove;
1520 
1521  // Need to check it here because we may not be
1522  // good after the read indicating that there
1523  // was nothing to read.
1524  if (!elementDataStrm.fail())
1525  {
1526  PointGeomSharedPtr vert(MemoryManager<PointGeom>::AllocateSharedPtr(m_meshDimension, edgeindx, xval, yval, zval));
1527 
1528  curve->m_points.push_back(vert);
1529  }
1530 
1531  }
1532  }
1533  catch(...)
1534  {
1535  NEKERROR(ErrorUtil::efatal,
1536  (std::string("Unable to read curve data for EDGE: ") + elementStr).c_str());
1537 
1538  }
1539 
1540  ASSERTL0(curve->m_points.size() == numPts,
1541  "Number of points specificed by attribute "
1542  "NUMPOINTS is different from number of points "
1543  "in list (edgeid = " +
1544  boost::lexical_cast<string>(edgeid));
1545 
1546  m_curvedEdges[edgeid] = curve;
1547 
1548  edgelement = edgelement->NextSiblingElement("E");
1549 
1550  } // end if-loop
1551 
1552  } // end while-loop
1553 
1554  TiXmlElement *facelement = field->FirstChildElement("F");
1555  int faceindx, faceid;
1556 
1557  while(facelement)
1558  {
1559  std::string face(facelement->ValueStr());
1560  ASSERTL0(face == "F", (std::string("Unknown 3D curve type: ") + face).c_str());
1561 
1562  /// Read id attribute.
1563  err = facelement->QueryIntAttribute("ID", &faceindx);
1564  ASSERTL0(err == TIXML_SUCCESS, "Unable to read curve attribute ID.");
1565 
1566  /// Read face id attribute.
1567  err = facelement->QueryIntAttribute("FACEID", &faceid);
1568  ASSERTL0(err == TIXML_SUCCESS, "Unable to read curve attribute FACEID.");
1569 
1570  /// Read text face element description.
1571  std::string elementStr;
1572  TiXmlNode* elementChild = facelement->FirstChild();
1573 
1574  while(elementChild)
1575  {
1576  // Accumulate all non-comment element data
1577  if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
1578  {
1579  elementStr += elementChild->ToText()->ValueStr();
1580  elementStr += " ";
1581  }
1582  elementChild = elementChild->NextSibling();
1583  }
1584 
1585  ASSERTL0(!elementStr.empty(), "Unable to read curve description body.");
1586 
1587  /// Parse out the element components corresponding to type of element.
1588  if(face == "F")
1589  {
1590  std::string typeStr = facelement->Attribute("TYPE");
1591  ASSERTL0(!typeStr.empty(), "TYPE must be specified in " "points definition");
1592  LibUtilities::PointsType type;
1593  const std::string* begStr = LibUtilities::kPointsTypeStr;
1594  const std::string* endStr = LibUtilities::kPointsTypeStr + LibUtilities::SIZE_PointsType;
1595  const std::string* ptsStr = std::find(begStr, endStr, typeStr);
1596 
1597  ASSERTL0(ptsStr != endStr, "Invalid points type.");
1598  type = (LibUtilities::PointsType)(ptsStr - begStr);
1599 
1600  std::string numptsStr = facelement->Attribute("NUMPOINTS");
1601  ASSERTL0(!numptsStr.empty(), "NUMPOINTS must be specified in points definition");
1602  int numPts=0;
1603  std::stringstream s;
1604  s << numptsStr;
1605  s >> numPts;
1606 
1607  CurveSharedPtr curve(MemoryManager<Curve>::AllocateSharedPtr(faceid, type));
1608 
1609  ASSERTL0(numPts >= 3, "NUMPOINTS for face must be greater than 2");
1610 
1611  if(numPts == 3)
1612  {
1613  ASSERTL0(ptsStr != endStr, "Invalid points type.");
1614  }
1615 
1616  // Read points (x, y, z)
1617  NekDouble xval, yval, zval;
1618  std::istringstream elementDataStrm(elementStr.c_str());
1619  try
1620  {
1621  while(!elementDataStrm.fail())
1622  {
1623  elementDataStrm >> xval >> yval >> zval;
1624 
1625  // Need to check it here because we
1626  // may not be good after the read
1627  // indicating that there was nothing
1628  // to read.
1629  if (!elementDataStrm.fail())
1630  {
1631  PointGeomSharedPtr vert(MemoryManager<PointGeom>::AllocateSharedPtr(m_meshDimension, faceindx, xval, yval, zval));
1632  curve->m_points.push_back(vert);
1633  }
1634  }
1635  }
1636  catch(...)
1637  {
1638  NEKERROR(ErrorUtil::efatal,
1639  (std::string("Unable to read curve data for FACE: ")
1640  + elementStr).c_str());
1641  }
1642  m_curvedFaces[faceid] = curve;
1643 
1644  facelement = facelement->NextSiblingElement("F");
1645 
1646  } // end if-loop
1647  } // end while-loop
1648  } // end of compressed else
1649  } // end of ReadCurves()
1650 
1651  /**
1652  *
1653  */
1654  void MeshGraph::ReadCurves(std::string &infilename)
1655  {
1656  TiXmlDocument doc(infilename);
1657  bool loadOkay = doc.LoadFile();
1658 
1659  std::stringstream errstr;
1660  errstr << "Unable to load file: " << infilename << std::endl;
1661  errstr << "Reason: " << doc.ErrorDesc() << std::endl;
1662  errstr << "Position: Line " << doc.ErrorRow() << ", Column " << doc.ErrorCol() << std::endl;
1663  ASSERTL0(loadOkay, errstr.str());
1664 
1665  ReadCurves(doc);
1666  }
1667 
1668  /**
1669  * @brief Populate a TinyXML document with a GEOMETRY tag inside the
1670  * NEKTAR tag.
1671  *
1672  * This routine will create a GEOMETRY XML tag which represents the
1673  * MeshGraph object. If a NEKTAR tag does not already exist, it will
1674  * create one. If a GEOMETRY block already exists inside the NEKTAR tag,
1675  * it will overwrite it.
1676  */
1677  void MeshGraph::WriteGeometry(TiXmlDocument &doc)
1678  {
1679  TiXmlElement *root = doc.FirstChildElement("NEKTAR");
1680  TiXmlElement *geomTag;
1681 
1682  // Try to find existing NEKTAR tag.
1683  if (!root)
1684  {
1685  root = new TiXmlElement("NEKTAR");
1686  doc.LinkEndChild(root);
1687 
1688  geomTag = new TiXmlElement("GEOMETRY");
1689  root->LinkEndChild(geomTag);
1690  }
1691  else
1692  {
1693  // Try to find existing GEOMETRY tag.
1694  geomTag = root->FirstChildElement("GEOMETRY");
1695 
1696  if (!geomTag)
1697  {
1698  geomTag = new TiXmlElement("GEOMETRY");
1699  root->LinkEndChild(geomTag);
1700  }
1701  }
1702 
1703  // Update attributes with dimensions.
1704  geomTag->SetAttribute("DIM", m_meshDimension);
1705  geomTag->SetAttribute("SPACE", m_spaceDimension);
1706 
1707  // Clear existing elements.
1708  geomTag->Clear();
1709 
1710  // Construct <VERTEX> block
1711  TiXmlElement *vertTag = new TiXmlElement("VERTEX");
1712  PointGeomMap::iterator pIt;
1713 
1714  for (pIt = m_vertSet.begin(); pIt != m_vertSet.end(); ++pIt)
1715  {
1716  stringstream s;
1717  s << scientific << setprecision(8)
1718  << (*pIt->second)(0) << " " << (*pIt->second)(1) << " "
1719  << (*pIt->second)(2);
1720  TiXmlElement * v = new TiXmlElement("V");
1721  v->SetAttribute("ID", pIt->second->GetVid());
1722  v->LinkEndChild(new TiXmlText(s.str()));
1723  vertTag->LinkEndChild(v);
1724  }
1725 
1726  geomTag->LinkEndChild(vertTag);
1727 
1728  // Construct <EDGE> or <ELEMENT> block
1729  TiXmlElement *edgeTag = new TiXmlElement(
1730  m_meshDimension == 1 ? "ELEMENT" : "EDGE");
1731  SegGeomMap::iterator sIt;
1732  string tag = m_meshDimension == 1 ? "S" : "E";
1733 
1734  for (sIt = m_segGeoms.begin(); sIt != m_segGeoms.end(); ++sIt)
1735  {
1736  stringstream s;
1737  SegGeomSharedPtr seg = sIt->second;
1738  s << seg->GetVid(0) << " " << seg->GetVid(1);
1739  TiXmlElement *e = new TiXmlElement(tag);
1740  e->SetAttribute("ID", sIt->first);
1741  e->LinkEndChild(new TiXmlText(s.str()));
1742  edgeTag->LinkEndChild(e);
1743  }
1744 
1745  geomTag->LinkEndChild(edgeTag);
1746 
1747  // Construct <FACE> or <ELEMENT> block
1748  if (m_meshDimension > 1)
1749  {
1750  TiXmlElement *faceTag = new TiXmlElement(
1751  m_meshDimension == 2 ? "ELEMENT" : "FACE");
1752 
1753  TriGeomMap::iterator tIt;
1754  tag = "T";
1755 
1756  for (tIt = m_triGeoms.begin(); tIt != m_triGeoms.end(); ++tIt)
1757  {
1758  stringstream s;
1759  TriGeomSharedPtr tri = tIt->second;
1760  s << tri->GetEid(0) << " " << tri->GetEid(1) << " "
1761  << tri->GetEid(2);
1762  TiXmlElement *t = new TiXmlElement(tag);
1763  t->SetAttribute("ID", tIt->first);
1764  t->LinkEndChild(new TiXmlText(s.str()));
1765  faceTag->LinkEndChild(t);
1766  }
1767 
1768  QuadGeomMap::iterator qIt;
1769  tag = "Q";
1770 
1771  for (qIt = m_quadGeoms.begin(); qIt != m_quadGeoms.end(); ++qIt)
1772  {
1773  stringstream s;
1774  QuadGeomSharedPtr quad = qIt->second;
1775  s << quad->GetEid(0) << " " << quad->GetEid(1) << " "
1776  << quad->GetEid(2) << " " << quad->GetEid(3);
1777  TiXmlElement *q = new TiXmlElement(tag);
1778  q->SetAttribute("ID", qIt->first);
1779  q->LinkEndChild(new TiXmlText(s.str()));
1780  faceTag->LinkEndChild(q);
1781  }
1782 
1783  geomTag->LinkEndChild(faceTag);
1784  }
1785 
1786  if (m_meshDimension > 2)
1787  {
1788  TiXmlElement *elmtTag = new TiXmlElement("ELEMENT");
1789 
1790  HexGeomMap::iterator hIt;
1791  tag = "H";
1792 
1793  for (hIt = m_hexGeoms.begin(); hIt != m_hexGeoms.end(); ++hIt)
1794  {
1795  stringstream s;
1796  HexGeomSharedPtr hex = hIt->second;
1797  s << hex->GetFid(0) << " " << hex->GetFid(1) << " "
1798  << hex->GetFid(2) << " " << hex->GetFid(3) << " "
1799  << hex->GetFid(4) << " " << hex->GetFid(5) << " ";
1800  TiXmlElement *h = new TiXmlElement(tag);
1801  h->SetAttribute("ID", hIt->first);
1802  h->LinkEndChild(new TiXmlText(s.str()));
1803  elmtTag->LinkEndChild(h);
1804  }
1805 
1806  PrismGeomMap::iterator rIt;
1807  tag = "R";
1808 
1809  for (rIt = m_prismGeoms.begin(); rIt != m_prismGeoms.end(); ++rIt)
1810  {
1811  stringstream s;
1812  PrismGeomSharedPtr prism = rIt->second;
1813  s << prism->GetFid(0) << " " << prism->GetFid(1) << " "
1814  << prism->GetFid(2) << " " << prism->GetFid(3) << " "
1815  << prism->GetFid(4) << " ";
1816  TiXmlElement *p = new TiXmlElement(tag);
1817  p->SetAttribute("ID", rIt->first);
1818  p->LinkEndChild(new TiXmlText(s.str()));
1819  elmtTag->LinkEndChild(p);
1820  }
1821 
1822  PyrGeomMap::iterator pIt;
1823  tag = "P";
1824 
1825  for (pIt = m_pyrGeoms.begin(); pIt != m_pyrGeoms.end(); ++pIt)
1826  {
1827  stringstream s;
1828  PyrGeomSharedPtr pyr = pIt->second;
1829  s << pyr->GetFid(0) << " " << pyr->GetFid(1) << " "
1830  << pyr->GetFid(2) << " " << pyr->GetFid(3) << " "
1831  << pyr->GetFid(4) << " ";
1832  TiXmlElement *p = new TiXmlElement(tag);
1833  p->SetAttribute("ID", pIt->first);
1834  p->LinkEndChild(new TiXmlText(s.str()));
1835  elmtTag->LinkEndChild(p);
1836  }
1837 
1838  TetGeomMap::iterator tIt;
1839  tag = "A";
1840 
1841  for (tIt = m_tetGeoms.begin(); tIt != m_tetGeoms.end(); ++tIt)
1842  {
1843  stringstream s;
1844  TetGeomSharedPtr tet = tIt->second;
1845  s << tet->GetFid(0) << " " << tet->GetFid(1) << " "
1846  << tet->GetFid(2) << " " << tet->GetFid(3) << " ";
1847  TiXmlElement *t = new TiXmlElement(tag);
1848  t->SetAttribute("ID", tIt->first);
1849  t->LinkEndChild(new TiXmlText(s.str()));
1850  elmtTag->LinkEndChild(t);
1851  }
1852 
1853  geomTag->LinkEndChild(elmtTag);
1854  }
1855 
1856  // Construct <CURVED> block
1857  TiXmlElement *curveTag = new TiXmlElement("CURVED");
1858  CurveMap::iterator curveIt;
1859  int curveId = 0;
1860 
1861  for (curveIt = m_curvedEdges.begin();
1862  curveIt != m_curvedEdges.end(); ++curveIt)
1863  {
1864  CurveSharedPtr curve = curveIt->second;
1865  TiXmlElement *c = new TiXmlElement("E");
1866  stringstream s;
1867  s.precision(8);
1868 
1869  for (int j = 0; j < curve->m_points.size(); ++j)
1870  {
1871  SpatialDomains::PointGeomSharedPtr p = curve->m_points[j];
1872  s << scientific << (*p)(0) << " " << (*p)(1) << " " << (*p)(2) << " ";
1873  }
1874 
1875  c->SetAttribute("ID", curveId++);
1876  c->SetAttribute("EDGEID", curve->m_curveID);
1877  c->SetAttribute("NUMPOINTS", curve->m_points.size());
1878  c->SetAttribute("TYPE", LibUtilities::kPointsTypeStr[curve->m_ptype]);
1879  c->LinkEndChild(new TiXmlText(s.str()));
1880  curveTag->LinkEndChild(c);
1881  }
1882 
1883  for (curveIt = m_curvedFaces.begin();
1884  curveIt != m_curvedFaces.end(); ++curveIt)
1885  {
1886  CurveSharedPtr curve = curveIt->second;
1887  TiXmlElement *c = new TiXmlElement("F");
1888  stringstream s;
1889  s.precision(8);
1890 
1891  for (int j = 0; j < curve->m_points.size(); ++j)
1892  {
1893  SpatialDomains::PointGeomSharedPtr p = curve->m_points[j];
1894  s << scientific << (*p)(0) << " " << (*p)(1) << " " << (*p)(2) << " ";
1895  }
1896 
1897  c->SetAttribute("ID", curveId++);
1898  c->SetAttribute("FACEID", curve->m_curveID);
1899  c->SetAttribute("NUMPOINTS", curve->m_points.size());
1900  c->SetAttribute("TYPE", LibUtilities::kPointsTypeStr[curve->m_ptype]);
1901  c->LinkEndChild(new TiXmlText(s.str()));
1902  curveTag->LinkEndChild(c);
1903  }
1904 
1905  geomTag->LinkEndChild(curveTag);
1906 
1907  // Construct <COMPOSITE> blocks
1908  TiXmlElement *compTag = new TiXmlElement("COMPOSITE");
1909  CompositeMap::iterator cIt;
1910 
1911  // Create a map that gets around the issue of mapping faces -> F and
1912  // edges -> E inside the tag.
1913  map<LibUtilities::ShapeType, pair<string, string> > compMap;
1914  compMap[LibUtilities::eSegment] = make_pair("S", "E");
1915  compMap[LibUtilities::eQuadrilateral] = make_pair("Q", "F");
1916  compMap[LibUtilities::eTriangle] = make_pair("T", "F");
1917  compMap[LibUtilities::eTetrahedron] = make_pair("A", "A");
1918  compMap[LibUtilities::ePyramid] = make_pair("P", "P");
1919  compMap[LibUtilities::ePrism] = make_pair("R", "R");
1920  compMap[LibUtilities::eHexahedron] = make_pair("H", "H");
1921 
1922  std::vector<unsigned int> idxList;
1923 
1924  for (cIt = m_meshComposites.begin(); cIt != m_meshComposites.end(); ++cIt)
1925  {
1926  stringstream s;
1927  TiXmlElement *c = new TiXmlElement("C");
1928 
1929  if (cIt->second->size() == 0)
1930  {
1931  continue;
1932  }
1933 
1934  GeometrySharedPtr firstGeom = cIt->second->at(0);
1935  int shapeDim = firstGeom->GetShapeDim();
1936  string tag = (shapeDim < m_meshDimension) ?
1937  compMap[firstGeom->GetShapeType()].second :
1938  compMap[firstGeom->GetShapeType()].first;
1939 
1940  idxList.clear();
1941  s << " " << tag << "[";
1942 
1943  for (int i = 0; i < cIt->second->size(); ++i)
1944  {
1945  idxList.push_back((*cIt->second)[i]->GetGlobalID());
1946  }
1947 
1948  s << ParseUtils::GenerateSeqString(idxList) << "] ";
1949 
1950  c->SetAttribute("ID", cIt->first);
1951  c->LinkEndChild(new TiXmlText(s.str()));
1952  compTag->LinkEndChild(c);
1953  }
1954 
1955  geomTag->LinkEndChild(compTag);
1956 
1957  // Construct <DOMAIN> block
1958  TiXmlElement *domTag = new TiXmlElement("DOMAIN");
1959  stringstream domString;
1960 
1961  // @todo Fix this to accomodate multi domain output
1962  idxList.clear();
1963  for (cIt = m_domain[0].begin(); cIt != m_domain[0].end(); ++cIt)
1964  {
1965  idxList.push_back(cIt->first);
1966  }
1967 
1968  domString << " C[" << ParseUtils::GenerateSeqString(idxList) << "] ";
1969  domTag->LinkEndChild(new TiXmlText(domString.str()));
1970  geomTag->LinkEndChild(domTag);
1971  }
1972 
1973  /**
1974  * @brief Write out an XML file containing the GEOMETRY block
1975  * representing this MeshGraph instance inside a NEKTAR tag.
1976  */
1977  void MeshGraph::WriteGeometry(std::string &outfilename)
1978  {
1979  // Create empty TinyXML document.
1980  TiXmlDocument doc;
1981  TiXmlDeclaration* decl = new TiXmlDeclaration( "1.0", "utf-8", "");
1982  doc.LinkEndChild(decl);
1983 
1984  // Write out geometry information.
1985  WriteGeometry(doc);
1986 
1987  // Save file.
1988  doc.SaveFile(outfilename);
1989  }
1990 
1991  void MeshGraph::SetDomainRange (
1992  NekDouble xmin, NekDouble xmax, NekDouble ymin,
1993  NekDouble ymax, NekDouble zmin, NekDouble zmax)
1994  {
1995  m_domainRange->m_checkShape = false;
1996 
1997  if(m_domainRange == NullDomainRangeShPtr)
1998  {
1999  m_domainRange = MemoryManager<DomainRange>::AllocateSharedPtr();
2000  m_domainRange->m_doXrange = true;
2001  }
2002 
2003  m_domainRange->m_xmin = xmin;
2004  m_domainRange->m_xmax = xmax;
2005 
2006  if(ymin == NekConstants::kNekUnsetDouble)
2007  {
2008  m_domainRange->m_doYrange = false;
2009  }
2010  else
2011  {
2012  m_domainRange->m_doYrange = true;
2013  m_domainRange->m_ymin = ymin;
2014  m_domainRange->m_ymax = ymax;
2015  }
2016 
2017  if(zmin == NekConstants::kNekUnsetDouble)
2018  {
2019  m_domainRange->m_doZrange = false;
2020  }
2021  else
2022  {
2023  m_domainRange->m_doZrange = true;
2024  m_domainRange->m_zmin = zmin;
2025  m_domainRange->m_zmax = zmax;
2026  }
2027  }
2028 
2029  bool MeshGraph::CheckRange(Geometry2D &geom)
2030  {
2031  bool returnval = true;
2032 
2033  if(m_domainRange != NullDomainRangeShPtr)
2034  {
2035  int nverts = geom.GetNumVerts();
2036  int coordim = geom.GetCoordim();
2037 
2038  // exclude elements outside x range if all vertices not in region
2039  if(m_domainRange->m_doXrange)
2040  {
2041  int ncnt_low = 0;
2042  int ncnt_up = 0;
2043  for(int i = 0; i < nverts; ++i)
2044  {
2045  NekDouble xval = (*geom.GetVertex(i))[0];
2046  if(xval < m_domainRange->m_xmin)
2047  {
2048  ncnt_low++;
2049  }
2050 
2051  if(xval > m_domainRange->m_xmax)
2052  {
2053  ncnt_up++;
2054  }
2055  }
2056 
2057  // check for all verts to be less or greater than
2058  // range so that if element spans thin range then
2059  // it is still included
2060  if((ncnt_up == nverts)||(ncnt_low == nverts))
2061  {
2062  returnval = false;
2063  }
2064  }
2065 
2066  // exclude elements outside y range if all vertices not in region
2067  if(m_domainRange->m_doYrange)
2068  {
2069  int ncnt_low = 0;
2070  int ncnt_up = 0;
2071  for(int i = 0; i < nverts; ++i)
2072  {
2073  NekDouble yval = (*geom.GetVertex(i))[1];
2074  if(yval < m_domainRange->m_ymin)
2075  {
2076  ncnt_low++;
2077  }
2078 
2079  if(yval > m_domainRange->m_ymax)
2080  {
2081  ncnt_up++;
2082  }
2083  }
2084 
2085  // check for all verts to be less or greater than
2086  // range so that if element spans thin range then
2087  // it is still included
2088  if((ncnt_up == nverts)||(ncnt_low == nverts))
2089  {
2090  returnval = false;
2091  }
2092  }
2093 
2094  if(coordim > 2)
2095  {
2096  // exclude elements outside z range if all vertices not in region
2097  if(m_domainRange->m_doZrange)
2098  {
2099  int ncnt_low = 0;
2100  int ncnt_up = 0;
2101 
2102  for(int i = 0; i < nverts; ++i)
2103  {
2104  NekDouble zval = (*geom.GetVertex(i))[2];
2105 
2106  if(zval < m_domainRange->m_zmin)
2107  {
2108  ncnt_low++;
2109  }
2110 
2111  if(zval > m_domainRange->m_zmax)
2112  {
2113  ncnt_up++;
2114  }
2115  }
2116 
2117  // check for all verts to be less or greater than
2118  // range so that if element spans thin range then
2119  // it is still included
2120  if((ncnt_up == nverts)||(ncnt_low == nverts))
2121  {
2122  returnval = false;
2123  }
2124  }
2125  }
2126  }
2127  return returnval;
2128  }
2129 
2130 
2131  /* Domain checker for 3D geometries */
2132  bool MeshGraph::CheckRange(Geometry3D &geom)
2133  {
2134  bool returnval = true;
2135 
2136  if(m_domainRange != NullDomainRangeShPtr)
2137  {
2138  int nverts = geom.GetNumVerts();
2139 
2140  if(m_domainRange->m_doXrange)
2141  {
2142  int ncnt_low = 0;
2143  int ncnt_up = 0;
2144 
2145  for(int i = 0; i < nverts; ++i)
2146  {
2147  NekDouble xval = (*geom.GetVertex(i))[0];
2148  if(xval < m_domainRange->m_xmin)
2149  {
2150  ncnt_low++;
2151  }
2152 
2153  if(xval > m_domainRange->m_xmax)
2154  {
2155  ncnt_up++;
2156  }
2157  }
2158 
2159  // check for all verts to be less or greater than
2160  // range so that if element spans thin range then
2161  // it is still included
2162  if((ncnt_up == nverts)||(ncnt_low == nverts))
2163  {
2164  returnval = false;
2165  }
2166  }
2167 
2168  if(m_domainRange->m_doYrange)
2169  {
2170  int ncnt_low = 0;
2171  int ncnt_up = 0;
2172  for(int i = 0; i < nverts; ++i)
2173  {
2174  NekDouble yval = (*geom.GetVertex(i))[1];
2175  if(yval < m_domainRange->m_ymin)
2176  {
2177  ncnt_low++;
2178  }
2179 
2180  if(yval > m_domainRange->m_ymax)
2181  {
2182  ncnt_up++;
2183  }
2184  }
2185 
2186  // check for all verts to be less or greater than
2187  // range so that if element spans thin range then
2188  // it is still included
2189  if((ncnt_up == nverts)||(ncnt_low == nverts))
2190  {
2191  returnval = false;
2192  }
2193  }
2194 
2195  if(m_domainRange->m_doZrange)
2196  {
2197  int ncnt_low = 0;
2198  int ncnt_up = 0;
2199  for(int i = 0; i < nverts; ++i)
2200  {
2201  NekDouble zval = (*geom.GetVertex(i))[2];
2202 
2203  if(zval < m_domainRange->m_zmin)
2204  {
2205  ncnt_low++;
2206  }
2207 
2208  if(zval > m_domainRange->m_zmax)
2209  {
2210  ncnt_up++;
2211  }
2212  }
2213 
2214  // check for all verts to be less or greater than
2215  // range so that if element spans thin range then
2216  // it is still included
2217  if((ncnt_up == nverts)||(ncnt_low == nverts))
2218  {
2219  returnval = false;
2220  }
2221  }
2222 
2223  if(m_domainRange->m_checkShape)
2224  {
2225  if(geom.GetShapeType() != m_domainRange->m_shapeType)
2226  {
2227  returnval = false;
2228  }
2229  }
2230 
2231  }
2232 
2233  return returnval;
2234  }
2235 
2236  /**
2237  *
2238  */
2239  GeometrySharedPtr MeshGraph::GetCompositeItem(int whichComposite, int whichItem)
2240  {
2241  GeometrySharedPtr returnval;
2242  bool error = false;
2243 
2244  if (whichComposite >= 0 && whichComposite < int(m_meshComposites.size()))
2245  {
2246  if (whichItem >= 0 && whichItem < int(m_meshComposites[whichComposite]->size()))
2247  {
2248  returnval = m_meshComposites[whichComposite]->at(whichItem);
2249  }
2250  else
2251  {
2252  error = true;
2253  }
2254  }
2255  else
2256  {
2257  error = true;
2258  }
2259 
2260  if (error)
2261  {
2262  std::ostringstream errStream;
2263  errStream << "Unable to access composite item [" << whichComposite << "][" << whichItem << "].";
2264 
2265  std::string testStr = errStream.str();
2266 
2267  NEKERROR(ErrorUtil::efatal, testStr.c_str());
2268  }
2269 
2270  return returnval;
2271  }
2272 
2273 
2274  /**
2275  *
2276  */
2277  void MeshGraph::GetCompositeList(const std::string &compositeStr, CompositeMap &compositeVector) const
2278  {
2279  // Parse the composites into a list.
2280  typedef vector<unsigned int> SeqVector;
2281  SeqVector seqVector;
2282  bool parseGood = ParseUtils::GenerateSeqVector(compositeStr.c_str(), seqVector);
2283 
2284  ASSERTL0(parseGood && !seqVector.empty(), (std::string("Unable to read composite index range: ") + compositeStr).c_str());
2285 
2286  SeqVector addedVector; // Vector of those composites already added to compositeVector;
2287  for (SeqVector::iterator iter = seqVector.begin(); iter != seqVector.end(); ++iter)
2288  {
2289  // Only add a new one if it does not already exist in vector.
2290  // Can't go back and delete with a vector, so prevent it from
2291  // being added in the first place.
2292  if (std::find(addedVector.begin(), addedVector.end(), *iter) == addedVector.end())
2293  {
2294 
2295  // If the composite listed is not found and we are working
2296  // on a partitioned mesh, silently ignore it.
2297  if (m_meshComposites.find(*iter) == m_meshComposites.end()
2298  && m_meshPartitioned)
2299  {
2300  continue;
2301  }
2302 
2303  addedVector.push_back(*iter);
2304  Composite composite = GetComposite(*iter);
2305  CompositeMap::iterator compIter;
2306  if (composite)
2307  {
2308  compositeVector[*iter] = composite;
2309  }
2310  else
2311  {
2312  char str[64];
2313  ::sprintf(str, "%d", *iter);
2314  NEKERROR(ErrorUtil::ewarning, (std::string("Undefined composite: ") + str).c_str());
2315 
2316  }
2317  }
2318  }
2319  }
2320 
2321 
2322  /**
2323  *
2324  */
2325  const ExpansionMap &MeshGraph::GetExpansions(const std::string variable)
2326  {
2327  ExpansionMapShPtr returnval;
2328 
2329  if(m_expansionMapShPtrMap.count(variable))
2330  {
2331  returnval = m_expansionMapShPtrMap.find(variable)->second;
2332  }
2333  else
2334  {
2335  if(m_expansionMapShPtrMap.count("DefaultVar") == 0)
2336  {
2337  NEKERROR(ErrorUtil::efatal, (std::string("Unable to find expansion vector definition for field: ")+variable).c_str());
2338  }
2339  returnval = m_expansionMapShPtrMap.find("DefaultVar")->second;
2340  m_expansionMapShPtrMap[variable] = returnval;
2341 
2342  NEKERROR(ErrorUtil::ewarning, (std::string("Using Default variable expansion definition for field: ")+variable).c_str());
2343  }
2344 
2345  return *returnval;
2346  }
2347 
2348 
2349  /**
2350  *
2351  */
2352  ExpansionShPtr MeshGraph::GetExpansion(GeometrySharedPtr geom, const std::string variable)
2353  {
2354  ExpansionMapIter iter;
2355  ExpansionMapShPtr expansionMap = m_expansionMapShPtrMap.find(variable)->second;
2356 
2357  iter = expansionMap->find(geom->GetGlobalID());
2358  ASSERTL1(iter != expansionMap->end(),
2359  "Could not find expansion " +
2360  boost::lexical_cast<string>(geom->GetGlobalID()) +
2361  " in expansion for variable " + variable);
2362  return iter->second;
2363  }
2364 
2365 
2366  /**
2367  *
2368  */
2369  void MeshGraph::SetExpansions(
2370  std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
2371  {
2372  int i, j, k, cnt, id;
2373  GeometrySharedPtr geom;
2374 
2375  ExpansionMapShPtr expansionMap;
2376 
2377  // Loop over fields and determine unique fields string and
2378  // declare whole expansion list
2379  for(i = 0; i < fielddef.size(); ++i)
2380  {
2381  for(j = 0; j < fielddef[i]->m_fields.size(); ++j)
2382  {
2383  std::string field = fielddef[i]->m_fields[j];
2384  if(m_expansionMapShPtrMap.count(field) == 0)
2385  {
2386  expansionMap = SetUpExpansionMap();
2387  m_expansionMapShPtrMap[field] = expansionMap;
2388 
2389  // check to see if DefaultVar also not set and
2390  // if so assign it to this expansion
2391  if(m_expansionMapShPtrMap.count("DefaultVar") == 0)
2392  {
2393  m_expansionMapShPtrMap["DefaultVar"] = expansionMap;
2394  }
2395  }
2396  }
2397  }
2398 
2399  // loop over all elements find the geometry shared ptr and
2400  // set up basiskey vector
2401  for(i = 0; i < fielddef.size(); ++i)
2402  {
2403  cnt = 0;
2404  std::vector<std::string> fields = fielddef[i]->m_fields;
2405  std::vector<unsigned int> nmodes = fielddef[i]->m_numModes;
2406  std::vector<LibUtilities::BasisType> basis = fielddef[i]->m_basis;
2407  bool pointDef = fielddef[i]->m_pointsDef;
2408  bool numPointDef = fielddef[i]->m_numPointsDef;
2409 
2410  // Check points and numpoints
2411  std::vector<unsigned int> npoints = fielddef[i]->m_numPoints;
2412  std::vector<LibUtilities::PointsType> points = fielddef[i]->m_points;
2413 
2414  bool UniOrder = fielddef[i]->m_uniOrder;
2415 
2416  for (j = 0; j < fielddef[i]->m_elementIDs.size(); ++j)
2417  {
2418 
2419  LibUtilities::BasisKeyVector bkeyvec;
2420  id = fielddef[i]->m_elementIDs[j];
2421 
2422  switch (fielddef[i]->m_shapeType)
2423  {
2424  case LibUtilities::eSegment:
2425  {
2426  if(m_segGeoms.count(fielddef[i]->m_elementIDs[j]) == 0)
2427  {
2428  // skip element likely from parallel read
2429  if(!UniOrder)
2430  {
2431  cnt++;
2432  cnt += fielddef[i]->m_numHomogeneousDir;
2433  }
2434  continue;
2435  }
2436  geom = m_segGeoms[fielddef[i]->m_elementIDs[j]];
2437 
2438  LibUtilities::PointsKey pkey(nmodes[cnt]+1, LibUtilities::eGaussLobattoLegendre);
2439 
2440  if(numPointDef&&pointDef)
2441  {
2442  const LibUtilities::PointsKey pkey1(npoints[cnt], points[0]);
2443  pkey = pkey1;
2444  }
2445  else if(!numPointDef&&pointDef)
2446  {
2447  const LibUtilities::PointsKey pkey1(nmodes[cnt]+1, points[0]);
2448  pkey = pkey1;
2449  }
2450  else if(numPointDef&&!pointDef)
2451  {
2452  const LibUtilities::PointsKey pkey1(npoints[cnt], LibUtilities::eGaussLobattoLegendre);
2453  pkey = pkey1;
2454  }
2455 
2456  LibUtilities::BasisKey bkey(basis[0], nmodes[cnt], pkey);
2457 
2458  if(!UniOrder)
2459  {
2460  cnt++;
2461  cnt += fielddef[i]->m_numHomogeneousDir;
2462  }
2463  bkeyvec.push_back(bkey);
2464  }
2465  break;
2466  case LibUtilities::eTriangle:
2467  {
2468  if(m_triGeoms.count(fielddef[i]->m_elementIDs[j]) == 0)
2469  {
2470  // skip element likely from parallel read
2471  if(!UniOrder)
2472  {
2473  cnt += 2;
2474  cnt += fielddef[i]->m_numHomogeneousDir;
2475  }
2476  continue;
2477  }
2478  geom = m_triGeoms[fielddef[i]->m_elementIDs[j]];
2479 
2480  LibUtilities::PointsKey pkey(nmodes[cnt]+1, LibUtilities::eGaussLobattoLegendre);
2481  if(numPointDef&&pointDef)
2482  {
2483  const LibUtilities::PointsKey pkey2(npoints[cnt], points[0]);
2484  pkey = pkey2;
2485  }
2486  else if(!numPointDef&&pointDef)
2487  {
2488  const LibUtilities::PointsKey pkey2(nmodes[cnt]+1, points[0]);
2489  pkey = pkey2;
2490  }
2491  else if(numPointDef&&!pointDef)
2492  {
2493  const LibUtilities::PointsKey pkey2(npoints[cnt], LibUtilities::eGaussLobattoLegendre);
2494  pkey = pkey2;
2495  }
2496  LibUtilities::BasisKey bkey(basis[0], nmodes[cnt], pkey);
2497 
2498  bkeyvec.push_back(bkey);
2499 
2500  LibUtilities::PointsKey pkey1(nmodes[cnt+1], LibUtilities::eGaussRadauMAlpha1Beta0);
2501  if(numPointDef&&pointDef)
2502  {
2503  const LibUtilities::PointsKey pkey2(npoints[cnt+1], points[1]);
2504  pkey1 = pkey2;
2505  }
2506  else if(!numPointDef&&pointDef)
2507  {
2508  const LibUtilities::PointsKey pkey2(nmodes[cnt+1], points[1]);
2509  pkey1 = pkey2;
2510  }
2511  else if(numPointDef&&!pointDef)
2512  {
2513  const LibUtilities::PointsKey pkey2(npoints[cnt+1], LibUtilities::eGaussRadauMAlpha1Beta0);
2514  pkey1 = pkey2;
2515  }
2516  LibUtilities::BasisKey bkey1(basis[1], nmodes[cnt+1], pkey1);
2517  bkeyvec.push_back(bkey1);
2518 
2519  if(!UniOrder)
2520  {
2521  cnt += 2;
2522  cnt += fielddef[i]->m_numHomogeneousDir;
2523  }
2524  }
2525  break;
2526  case LibUtilities::eQuadrilateral:
2527  {
2528  if(m_quadGeoms.count(fielddef[i]->m_elementIDs[j]) == 0)
2529  {
2530  // skip element likely from parallel read
2531  if(!UniOrder)
2532  {
2533  cnt += 2;
2534  cnt += fielddef[i]->m_numHomogeneousDir;
2535  }
2536  continue;
2537  }
2538 
2539  geom = m_quadGeoms[fielddef[i]->m_elementIDs[j]];
2540 
2541  for(int b = 0; b < 2; ++b)
2542  {
2543  LibUtilities::PointsKey pkey(nmodes[cnt+b]+1, LibUtilities::eGaussLobattoLegendre);
2544 
2545  if(numPointDef&&pointDef)
2546  {
2547  const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
2548  pkey = pkey2;
2549  }
2550  else if(!numPointDef&&pointDef)
2551  {
2552  const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
2553  pkey = pkey2;
2554  }
2555  else if(numPointDef&&!pointDef)
2556  {
2557  const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
2558  pkey = pkey2;
2559  }
2560  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2561  bkeyvec.push_back(bkey);
2562  }
2563 
2564  if(!UniOrder)
2565  {
2566  cnt += 2;
2567  cnt += fielddef[i]->m_numHomogeneousDir;
2568  }
2569  }
2570  break;
2571 
2572  case LibUtilities::eTetrahedron:
2573  {
2574  k = fielddef[i]->m_elementIDs[j];
2575 
2576  // allow for possibility that fielddef is
2577  // larger than m_graph which can happen in
2578  // parallel runs
2579  if(m_tetGeoms.count(k) == 0)
2580  {
2581  if(!UniOrder)
2582  {
2583  cnt += 3;
2584  }
2585  continue;
2586  }
2587  geom = m_tetGeoms[k];
2588 
2589  {
2590  LibUtilities::PointsKey pkey(nmodes[cnt]+1, LibUtilities::eGaussLobattoLegendre);
2591 
2592  if(numPointDef&&pointDef)
2593  {
2594  const LibUtilities::PointsKey pkey2(npoints[cnt],points[0]);
2595  pkey = pkey2;
2596  }
2597  else if(!numPointDef&&pointDef)
2598  {
2599  const LibUtilities::PointsKey pkey2(nmodes[cnt]+1,points[0]);
2600  pkey = pkey2;
2601  }
2602  else if(numPointDef&&!pointDef)
2603  {
2604  const LibUtilities::PointsKey pkey2(npoints[cnt],LibUtilities::eGaussLobattoLegendre);
2605  pkey = pkey2;
2606  }
2607 
2608  LibUtilities::BasisKey bkey(basis[0],nmodes[cnt],pkey);
2609 
2610  bkeyvec.push_back(bkey);
2611  }
2612  {
2613  LibUtilities::PointsKey pkey(nmodes[cnt+1], LibUtilities::eGaussRadauMAlpha1Beta0);
2614 
2615  if(numPointDef&&pointDef)
2616  {
2617  const LibUtilities::PointsKey pkey2(npoints[cnt+1],points[1]);
2618  pkey = pkey2;
2619  }
2620  else if(!numPointDef&&pointDef)
2621  {
2622  const LibUtilities::PointsKey pkey2(nmodes[cnt+1]+1,points[1]);
2623  pkey = pkey2;
2624  }
2625  else if(numPointDef&&!pointDef)
2626  {
2627  const LibUtilities::PointsKey pkey2(npoints[cnt+1],LibUtilities::eGaussRadauMAlpha1Beta0);
2628  pkey = pkey2;
2629  }
2630 
2631  LibUtilities::BasisKey bkey(basis[1],nmodes[cnt+1],pkey);
2632 
2633  bkeyvec.push_back(bkey);
2634  }
2635 
2636  {
2637  LibUtilities::PointsKey pkey(nmodes[cnt+2], LibUtilities::eGaussRadauMAlpha2Beta0);
2638 
2639  if(numPointDef&&pointDef)
2640  {
2641  const LibUtilities::PointsKey pkey2(npoints[cnt+2],points[2]);
2642  pkey = pkey2;
2643  }
2644  else if(!numPointDef&&pointDef)
2645  {
2646  const LibUtilities::PointsKey pkey2(nmodes[cnt+2]+1,points[2]);
2647  pkey = pkey2;
2648  }
2649  else if(numPointDef&&!pointDef)
2650  {
2651  const LibUtilities::PointsKey pkey2(npoints[cnt+2],LibUtilities::eGaussRadauMAlpha1Beta0);
2652  pkey = pkey2;
2653  }
2654 
2655  LibUtilities::BasisKey bkey(basis[2],nmodes[cnt+2],pkey);
2656 
2657  bkeyvec.push_back(bkey);
2658  }
2659 
2660  if(!UniOrder)
2661  {
2662  cnt += 3;
2663  }
2664  }
2665  break;
2666  case LibUtilities::ePrism:
2667  {
2668  k = fielddef[i]->m_elementIDs[j];
2669  if(m_prismGeoms.count(k) == 0)
2670  {
2671  if(!UniOrder)
2672  {
2673  cnt += 3;
2674  }
2675  continue;
2676  }
2677  geom = m_prismGeoms[k];
2678 
2679  for(int b = 0; b < 2; ++b)
2680  {
2681  LibUtilities::PointsKey pkey(nmodes[cnt+b]+1,LibUtilities::eGaussLobattoLegendre);
2682 
2683  if(numPointDef&&pointDef)
2684  {
2685  const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
2686  pkey = pkey2;
2687  }
2688  else if(!numPointDef&&pointDef)
2689  {
2690  const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
2691  pkey = pkey2;
2692  }
2693  else if(numPointDef&&!pointDef)
2694  {
2695  const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
2696  pkey = pkey2;
2697  }
2698 
2699  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2700  bkeyvec.push_back(bkey);
2701  }
2702 
2703  {
2704  LibUtilities::PointsKey pkey(nmodes[cnt+2],LibUtilities::eGaussRadauMAlpha1Beta0);
2705 
2706  if(numPointDef&&pointDef)
2707  {
2708  const LibUtilities::PointsKey pkey2(npoints[cnt+2],points[2]);
2709  pkey = pkey2;
2710  }
2711  else if(!numPointDef&&pointDef)
2712  {
2713  const LibUtilities::PointsKey pkey2(nmodes[cnt+2]+1,points[2]);
2714  pkey = pkey2;
2715  }
2716  else if(numPointDef&&!pointDef)
2717  {
2718  const LibUtilities::PointsKey pkey2(npoints[cnt+2],LibUtilities::eGaussLobattoLegendre);
2719  pkey = pkey2;
2720  }
2721 
2722  LibUtilities::BasisKey bkey(basis[2],nmodes[cnt+2],pkey);
2723  bkeyvec.push_back(bkey);
2724  }
2725 
2726  if(!UniOrder)
2727  {
2728  cnt += 3;
2729  }
2730  }
2731  break;
2732  case LibUtilities::ePyramid:
2733  {
2734  k = fielddef[i]->m_elementIDs[j];
2735  ASSERTL0(m_pyrGeoms.find(k) != m_pyrGeoms.end(),
2736  "Failed to find geometry with same global id");
2737  geom = m_pyrGeoms[k];
2738 
2739 
2740  for(int b = 0; b < 2; ++b)
2741  {
2742  LibUtilities::PointsKey pkey(nmodes[cnt+b]+1,LibUtilities::eGaussLobattoLegendre);
2743 
2744  if(numPointDef&&pointDef)
2745  {
2746  const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
2747  pkey = pkey2;
2748  }
2749  else if(!numPointDef&&pointDef)
2750  {
2751  const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
2752  pkey = pkey2;
2753  }
2754  else if(numPointDef&&!pointDef)
2755  {
2756  const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
2757  pkey = pkey2;
2758  }
2759 
2760  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2761  bkeyvec.push_back(bkey);
2762  }
2763 
2764  {
2765  LibUtilities::PointsKey pkey(nmodes[cnt+2],LibUtilities::eGaussRadauMAlpha2Beta0);
2766 
2767  if(numPointDef&&pointDef)
2768  {
2769  const LibUtilities::PointsKey pkey2(npoints[cnt+2],points[2]);
2770  pkey = pkey2;
2771  }
2772  else if(!numPointDef&&pointDef)
2773  {
2774  const LibUtilities::PointsKey pkey2(nmodes[cnt+2]+1,points[2]);
2775  pkey = pkey2;
2776  }
2777  else if(numPointDef&&!pointDef)
2778  {
2779  const LibUtilities::PointsKey pkey2(npoints[cnt+2],LibUtilities::eGaussLobattoLegendre);
2780  pkey = pkey2;
2781  }
2782 
2783  LibUtilities::BasisKey bkey(basis[2],nmodes[cnt+2],pkey);
2784  bkeyvec.push_back(bkey);
2785  }
2786 
2787  if(!UniOrder)
2788  {
2789  cnt += 3;
2790  }
2791  }
2792  break;
2793  case LibUtilities::eHexahedron:
2794  {
2795  k = fielddef[i]->m_elementIDs[j];
2796  if(m_hexGeoms.count(k) == 0)
2797  {
2798  if(!UniOrder)
2799  {
2800  cnt += 3;
2801  }
2802  continue;
2803  }
2804 
2805  geom = m_hexGeoms[k];
2806 
2807  for(int b = 0; b < 3; ++b)
2808  {
2809  LibUtilities::PointsKey pkey(nmodes[cnt+b],LibUtilities::eGaussLobattoLegendre);
2810 
2811  if(numPointDef&&pointDef)
2812  {
2813  const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
2814  pkey = pkey2;
2815  }
2816  else if(!numPointDef&&pointDef)
2817  {
2818  const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
2819  pkey = pkey2;
2820  }
2821  else if(numPointDef&&!pointDef)
2822  {
2823  const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
2824  pkey = pkey2;
2825  }
2826 
2827  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2828  bkeyvec.push_back(bkey);
2829  }
2830 
2831  if(!UniOrder)
2832  {
2833  cnt += 3;
2834  }
2835  }
2836  break;
2837  default:
2838  ASSERTL0(false,"Need to set up for pyramid and prism 3D Expansions");
2839  break;
2840  }
2841 
2842  for(k = 0; k < fields.size(); ++k)
2843  {
2844  expansionMap = m_expansionMapShPtrMap.find(fields[k])->second;
2845  if((*expansionMap).find(id) != (*expansionMap).end())
2846  {
2847  (*expansionMap)[id]->m_geomShPtr = geom;
2848  (*expansionMap)[id]->m_basisKeyVector = bkeyvec;
2849  }
2850  }
2851  }
2852  }
2853  }
2854 
2855 
2856  /**
2857  *
2858  */
2859  void MeshGraph::SetExpansions(
2860  std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef,
2861  std::vector< std::vector<LibUtilities::PointsType> > &pointstype)
2862  {
2863  int i,j,k,cnt,id;
2864  GeometrySharedPtr geom;
2865 
2866  ExpansionMapShPtr expansionMap;
2867 
2868  // Loop over fields and determine unique fields string and
2869  // declare whole expansion list
2870  for(i = 0; i < fielddef.size(); ++i)
2871  {
2872  for(j = 0; j < fielddef[i]->m_fields.size(); ++j)
2873  {
2874  std::string field = fielddef[i]->m_fields[j];
2875  if(m_expansionMapShPtrMap.count(field) == 0)
2876  {
2877  expansionMap = SetUpExpansionMap();
2878  m_expansionMapShPtrMap[field] = expansionMap;
2879 
2880  // check to see if DefaultVar also not set and
2881  // if so assign it to this expansion
2882  if(m_expansionMapShPtrMap.count("DefaultVar") == 0)
2883  {
2884  m_expansionMapShPtrMap["DefaultVar"] = expansionMap;
2885  }
2886  }
2887  }
2888  }
2889 
2890 
2891  // loop over all elements find the geometry shared ptr and
2892  // set up basiskey vector
2893  for(i = 0; i < fielddef.size(); ++i)
2894  {
2895  cnt = 0;
2896  std::vector<std::string> fields = fielddef[i]->m_fields;
2897  std::vector<unsigned int> nmodes = fielddef[i]->m_numModes;
2898  std::vector<LibUtilities::BasisType> basis = fielddef[i]->m_basis;
2899  bool UniOrder = fielddef[i]->m_uniOrder;
2900 
2901  for(j = 0; j < fielddef[i]->m_elementIDs.size(); ++j)
2902  {
2903  LibUtilities::BasisKeyVector bkeyvec;
2904  id = fielddef[i]->m_elementIDs[j];
2905 
2906  switch(fielddef[i]->m_shapeType)
2907  {
2908  case LibUtilities::eSegment:
2909  {
2910  k = fielddef[i]->m_elementIDs[j];
2911  ASSERTL0(m_segGeoms.find(k) != m_segGeoms.end(),
2912  "Failed to find geometry with same global id.");
2913  geom = m_segGeoms[k];
2914 
2915  const LibUtilities::PointsKey pkey(nmodes[cnt], pointstype[i][0]);
2916  LibUtilities::BasisKey bkey(basis[0], nmodes[cnt], pkey);
2917  if(!UniOrder)
2918  {
2919  cnt++;
2920  cnt += fielddef[i]->m_numHomogeneousDir;
2921  }
2922  bkeyvec.push_back(bkey);
2923  }
2924  break;
2925  case LibUtilities::eTriangle:
2926  {
2927  k = fielddef[i]->m_elementIDs[j];
2928  ASSERTL0(m_triGeoms.find(k) != m_triGeoms.end(),
2929  "Failed to find geometry with same global id.");
2930  geom = m_triGeoms[k];
2931  for(int b = 0; b < 2; ++b)
2932  {
2933  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
2934  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2935  bkeyvec.push_back(bkey);
2936  }
2937 
2938  if(!UniOrder)
2939  {
2940  cnt += 2;
2941  cnt += fielddef[i]->m_numHomogeneousDir;
2942  }
2943  }
2944  break;
2945  case LibUtilities::eQuadrilateral:
2946  {
2947  k = fielddef[i]->m_elementIDs[j];
2948  ASSERTL0(m_quadGeoms.find(k) != m_quadGeoms.end(),
2949  "Failed to find geometry with same global id");
2950  geom = m_quadGeoms[k];
2951 
2952  for(int b = 0; b < 2; ++b)
2953  {
2954  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
2955  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2956  bkeyvec.push_back(bkey);
2957  }
2958 
2959  if(!UniOrder)
2960  {
2961  cnt += 2;
2962  cnt += fielddef[i]->m_numHomogeneousDir;
2963  }
2964  }
2965  break;
2966  case LibUtilities::eTetrahedron:
2967  {
2968  k = fielddef[i]->m_elementIDs[j];
2969  ASSERTL0(m_tetGeoms.find(k) != m_tetGeoms.end(),
2970  "Failed to find geometry with same global id");
2971  geom = m_tetGeoms[k];
2972 
2973  for(int b = 0; b < 3; ++b)
2974  {
2975  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
2976  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2977  bkeyvec.push_back(bkey);
2978  }
2979 
2980  if(!UniOrder)
2981  {
2982  cnt += 3;
2983  }
2984  }
2985  break;
2986  case LibUtilities::ePyramid:
2987  {
2988  k = fielddef[i]->m_elementIDs[j];
2989  ASSERTL0(m_pyrGeoms.find(k) != m_pyrGeoms.end(),
2990  "Failed to find geometry with same global id");
2991  geom = m_pyrGeoms[k];
2992 
2993  for(int b = 0; b < 3; ++b)
2994  {
2995  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
2996  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2997  bkeyvec.push_back(bkey);
2998  }
2999 
3000  if(!UniOrder)
3001  {
3002  cnt += 3;
3003  }
3004  }
3005  break;
3006  case LibUtilities::ePrism:
3007  {
3008  k = fielddef[i]->m_elementIDs[j];
3009  ASSERTL0(m_prismGeoms.find(k) != m_prismGeoms.end(),
3010  "Failed to find geometry with same global id");
3011  geom = m_prismGeoms[k];
3012 
3013  for(int b = 0; b < 3; ++b)
3014  {
3015  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
3016  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
3017  bkeyvec.push_back(bkey);
3018  }
3019 
3020  if(!UniOrder)
3021  {
3022  cnt += 3;
3023  }
3024  }
3025  break;
3026  case LibUtilities::eHexahedron:
3027  {
3028  k = fielddef[i]->m_elementIDs[j];
3029  ASSERTL0(m_hexGeoms.find(k) != m_hexGeoms.end(),
3030  "Failed to find geometry with same global id");
3031  geom = m_hexGeoms[k];
3032 
3033  for(int b = 0; b < 3; ++b)
3034  {
3035  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
3036  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
3037  bkeyvec.push_back(bkey);
3038  }
3039 
3040  if(!UniOrder)
3041  {
3042  cnt += 3;
3043  }
3044  }
3045  break;
3046  default:
3047  ASSERTL0(false,"Need to set up for pyramid and prism 3D Expansions");
3048  break;
3049  }
3050 
3051  for(k = 0; k < fields.size(); ++k)
3052  {
3053  expansionMap = m_expansionMapShPtrMap.find(fields[k])->second;
3054  if((*expansionMap).find(id) != (*expansionMap).end())
3055  {
3056  (*expansionMap)[id]->m_geomShPtr = geom;
3057  (*expansionMap)[id]->m_basisKeyVector = bkeyvec;
3058  }
3059  }
3060  }
3061  }
3062  }
3063 
3064  /**
3065  * \brief Reset all points keys to have equispaced points with
3066  * optional arguemt of \a npoints which redefines how many
3067  * points are to be used.
3068  */
3069  void MeshGraph::SetExpansionsToEvenlySpacedPoints(int npoints)
3070  {
3071  ExpansionMapShPtrMapIter it;
3072 
3073  // iterate over all defined expansions
3074  for(it = m_expansionMapShPtrMap.begin(); it != m_expansionMapShPtrMap.end(); ++it)
3075  {
3076  ExpansionMapIter expIt;
3077 
3078  for(expIt = it->second->begin(); expIt != it->second->end(); ++expIt)
3079  {
3080  for(int i = 0; i < expIt->second->m_basisKeyVector.size(); ++i)
3081  {
3082  LibUtilities::BasisKey bkeyold = expIt->second->m_basisKeyVector[i];
3083 
3084  int npts;
3085 
3086  if(npoints) // use input
3087  {
3088  npts = npoints;
3089  }
3090  else
3091  {
3092  npts = bkeyold.GetNumModes();
3093  }
3094 
3095 
3096  const LibUtilities::PointsKey pkey(npts,LibUtilities::ePolyEvenlySpaced);
3097  LibUtilities::BasisKey bkeynew(bkeyold.GetBasisType(),bkeyold.GetNumModes(), pkey);
3098  expIt->second->m_basisKeyVector[i] = bkeynew;
3099 
3100  }
3101  }
3102  }
3103  }
3104 
3105  /**
3106  * \brief Reset all points keys to have expansion order of \a
3107  * nmodes. we keep the point distribution the same and make
3108  * the number of points the same difference from the number
3109  * of modes as the original expansion definition
3110  */
3111  void MeshGraph::SetExpansionsToPolyOrder(int nmodes)
3112  {
3113  ExpansionMapShPtrMapIter it;
3114 
3115  // iterate over all defined expansions
3116  for(it = m_expansionMapShPtrMap.begin(); it != m_expansionMapShPtrMap.end(); ++it)
3117  {
3118  ExpansionMapIter expIt;
3119 
3120  for(expIt = it->second->begin(); expIt != it->second->end(); ++expIt)
3121  {
3122  for(int i = 0; i < expIt->second->m_basisKeyVector.size(); ++i)
3123  {
3124  LibUtilities::BasisKey bkeyold = expIt->second->m_basisKeyVector[i];
3125 
3126  int npts = nmodes + (bkeyold.GetNumPoints() - bkeyold.GetNumModes());
3127 
3128  const LibUtilities::PointsKey pkey(npts,bkeyold.GetPointsType());
3129  LibUtilities::BasisKey bkeynew(bkeyold.GetBasisType(),nmodes, pkey);
3130  expIt->second->m_basisKeyVector[i] = bkeynew;
3131 
3132  }
3133  }
3134  }
3135  }
3136 
3137 
3138  /**
3139  * \brief Reset all points keys to have expansion order of \a
3140  * nmodes. we keep the point distribution the same and make
3141  * the number of points the same difference from the number
3142  * of modes as the original expansion definition
3143  */
3144  void MeshGraph::SetExpansionsToPointOrder(int npts)
3145  {
3146  ExpansionMapShPtrMapIter it;
3147 
3148  // iterate over all defined expansions
3149  for (it = m_expansionMapShPtrMap.begin();
3150  it != m_expansionMapShPtrMap.end();
3151  ++it)
3152  {
3153  ExpansionMapIter expIt;
3154 
3155  for (expIt = it->second->begin();
3156  expIt != it->second->end();
3157  ++expIt)
3158  {
3159  for(int i = 0;
3160  i < expIt->second->m_basisKeyVector.size();
3161  ++i)
3162  {
3163  LibUtilities::BasisKey bkeyold =
3164  expIt->second->m_basisKeyVector[i];
3165 
3166  const LibUtilities::PointsKey pkey(
3167  npts, bkeyold.GetPointsType());
3168 
3169  LibUtilities::BasisKey bkeynew(bkeyold.GetBasisType(),
3170  bkeyold.GetNumModes(),
3171  pkey);
3172  expIt->second->m_basisKeyVector[i] = bkeynew;
3173  }
3174  }
3175  }
3176  }
3177 
3178 
3179  /**
3180  * For each element of shape given by \a shape in field \a
3181  * var, replace the current BasisKeyVector describing the
3182  * expansion in each dimension, with the one provided by \a
3183  * keys.
3184  *
3185  * @TODO: Allow selection of elements through a CompositeVector,
3186  * as well as by type.
3187  *
3188  * @param shape The shape of elements to be changed.
3189  * @param keys The new basis vector to apply to those elements.
3190  */
3191  void MeshGraph::SetBasisKey(LibUtilities::ShapeType shape,
3192  LibUtilities::BasisKeyVector &keys,
3193  std::string var)
3194  {
3195  ExpansionMapIter elemIter;
3196 
3197  ExpansionMapShPtr expansionMap = m_expansionMapShPtrMap.find(var)->second;
3198 
3199  for (elemIter = expansionMap->begin(); elemIter != expansionMap->end(); ++elemIter)
3200  {
3201  if ((elemIter->second)->m_geomShPtr->GetShapeType() == shape)
3202  {
3203  (elemIter->second)->m_basisKeyVector = keys;
3204  }
3205  }
3206  }
3207 
3208 
3209  /**
3210  *
3211  */
3212  LibUtilities::BasisKeyVector MeshGraph::DefineBasisKeyFromExpansionType(
3213  GeometrySharedPtr in,
3214  ExpansionType type,
3215  const int nummodes)
3216  {
3217  LibUtilities::BasisKeyVector returnval;
3218 
3219  LibUtilities::ShapeType shape= in->GetShapeType();
3220 
3221  int quadoffset = 1;
3222  switch(type)
3223  {
3224  case eModified:
3225  case eModifiedGLLRadau10:
3226  quadoffset = 1;
3227  break;
3228  case eModifiedQuadPlus1:
3229  quadoffset = 2;
3230  break;
3231  case eModifiedQuadPlus2:
3232  quadoffset = 3;
3233  break;
3234  default:
3235  break;
3236  }
3237 
3238  switch(type)
3239  {
3240  case eModified:
3241  case eModifiedQuadPlus1:
3242  case eModifiedQuadPlus2:
3243  case eModifiedGLLRadau10:
3244  {
3245  switch (shape)
3246  {
3247  case LibUtilities::eSegment:
3248  {
3249  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3250  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3251  returnval.push_back(bkey);
3252  }
3253  break;
3254  case LibUtilities::eQuadrilateral:
3255  {
3256  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3257  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3258  returnval.push_back(bkey);
3259  returnval.push_back(bkey);
3260  }
3261  break;
3262  case LibUtilities::eHexahedron:
3263  {
3264  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3265  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3266  returnval.push_back(bkey);
3267  returnval.push_back(bkey);
3268  returnval.push_back(bkey);
3269  }
3270  break;
3271  case LibUtilities::eTriangle:
3272  {
3273  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3274  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3275  returnval.push_back(bkey);
3276 
3277  const LibUtilities::PointsKey pkey1(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha1Beta0);
3278  LibUtilities::BasisKey bkey1(LibUtilities::eModified_B, nummodes, pkey1);
3279 
3280  returnval.push_back(bkey1);
3281  }
3282  break;
3283  case LibUtilities::eTetrahedron:
3284  {
3285  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3286  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3287  returnval.push_back(bkey);
3288 
3289  const LibUtilities::PointsKey pkey1(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha1Beta0);
3290  LibUtilities::BasisKey bkey1(LibUtilities::eModified_B, nummodes, pkey1);
3291  returnval.push_back(bkey1);
3292 
3293  if(type == eModifiedGLLRadau10)
3294  {
3295  const LibUtilities::PointsKey pkey2(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha1Beta0);
3296  LibUtilities::BasisKey bkey2(LibUtilities::eModified_C, nummodes, pkey2);
3297  returnval.push_back(bkey2);
3298  }
3299  else
3300  {
3301  const LibUtilities::PointsKey pkey2(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha2Beta0);
3302  LibUtilities::BasisKey bkey2(LibUtilities::eModified_C, nummodes, pkey2);
3303  returnval.push_back(bkey2);
3304  }
3305  }
3306  break;
3307  case LibUtilities::ePyramid:
3308  {
3309  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3310  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3311  returnval.push_back(bkey);
3312  returnval.push_back(bkey);
3313 
3314  const LibUtilities::PointsKey pkey1(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha2Beta0);
3315  LibUtilities::BasisKey bkey1(LibUtilities::eModified_C, nummodes, pkey1);
3316  returnval.push_back(bkey1);
3317  }
3318  break;
3319  case LibUtilities::ePrism:
3320  {
3321  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3322  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3323  returnval.push_back(bkey);
3324  returnval.push_back(bkey);
3325 
3326  const LibUtilities::PointsKey pkey1(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha1Beta0);
3327  LibUtilities::BasisKey bkey1(LibUtilities::eModified_B, nummodes, pkey1);
3328  returnval.push_back(bkey1);
3329 
3330  }
3331  break;
3332  default:
3333  {
3334  ASSERTL0(false,"Expansion not defined in switch for this shape");
3335  }
3336  break;
3337  }
3338  }
3339  break;
3340 
3341  case eGLL_Lagrange:
3342  {
3343  switch(shape)
3344  {
3345  case LibUtilities::eSegment:
3346  {
3347  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3348  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3349  returnval.push_back(bkey);
3350  }
3351  break;
3352  case LibUtilities::eQuadrilateral:
3353  {
3354  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3355  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3356  returnval.push_back(bkey);
3357  returnval.push_back(bkey);
3358  }
3359  break;
3360  case LibUtilities::eTriangle: // define with corrects points key
3361  // and change to Ortho on construction
3362  {
3363  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3364  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3365  returnval.push_back(bkey);
3366 
3367  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussRadauMAlpha1Beta0);
3368  LibUtilities::BasisKey bkey1(LibUtilities::eOrtho_B, nummodes, pkey1);
3369  returnval.push_back(bkey1);
3370  }
3371  break;
3372  case LibUtilities::eHexahedron:
3373  {
3374  const LibUtilities::PointsKey pkey(nummodes+1,LibUtilities::eGaussLobattoLegendre);
3375  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3376 
3377  returnval.push_back(bkey);
3378  returnval.push_back(bkey);
3379  returnval.push_back(bkey);
3380  }
3381  break;
3382  default:
3383  {
3384  ASSERTL0(false, "Expansion not defined in switch for this shape");
3385  }
3386  break;
3387  }
3388  }
3389  break;
3390 
3391  case eGauss_Lagrange:
3392  {
3393  switch (shape)
3394  {
3395  case LibUtilities::eSegment:
3396  {
3397  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussLegendre);
3398  LibUtilities::BasisKey bkey(LibUtilities::eGauss_Lagrange, nummodes, pkey);
3399 
3400  returnval.push_back(bkey);
3401  }
3402  break;
3403  case LibUtilities::eQuadrilateral:
3404  {
3405  const LibUtilities::PointsKey pkey(nummodes,LibUtilities::eGaussGaussLegendre);
3406  LibUtilities::BasisKey bkey(LibUtilities::eGauss_Lagrange, nummodes, pkey);
3407 
3408  returnval.push_back(bkey);
3409  returnval.push_back(bkey);
3410  }
3411  break;
3412  case LibUtilities::eHexahedron:
3413  {
3414  const LibUtilities::PointsKey pkey(nummodes,LibUtilities::eGaussGaussLegendre);
3415  LibUtilities::BasisKey bkey(LibUtilities::eGauss_Lagrange, nummodes, pkey);
3416 
3417  returnval.push_back(bkey);
3418  returnval.push_back(bkey);
3419  returnval.push_back(bkey);
3420  }
3421  break;
3422  default:
3423  {
3424  ASSERTL0(false, "Expansion not defined in switch for this shape");
3425  }
3426  break;
3427  }
3428  }
3429  break;
3430 
3431  case eOrthogonal:
3432  {
3433  switch (shape)
3434  {
3435  case LibUtilities::eSegment:
3436  {
3437  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3438  LibUtilities::BasisKey bkey(LibUtilities::eOrtho_A, nummodes, pkey);
3439 
3440  returnval.push_back(bkey);
3441  }
3442  break;
3443  case LibUtilities::eTriangle:
3444  {
3445  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3446  LibUtilities::BasisKey bkey(LibUtilities::eOrtho_A, nummodes, pkey);
3447 
3448  returnval.push_back(bkey);
3449 
3450  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussRadauMAlpha1Beta0);
3451  LibUtilities::BasisKey bkey1(LibUtilities::eOrtho_B, nummodes, pkey1);
3452 
3453  returnval.push_back(bkey1);
3454  }
3455  break;
3456  case LibUtilities::eQuadrilateral:
3457  {
3458  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3459  LibUtilities::BasisKey bkey(LibUtilities::eOrtho_A, nummodes, pkey);
3460 
3461  returnval.push_back(bkey);
3462  returnval.push_back(bkey);
3463  }
3464  break;
3465  case LibUtilities::eTetrahedron:
3466  {
3467  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3468  LibUtilities::BasisKey bkey(LibUtilities::eOrtho_A, nummodes, pkey);
3469 
3470  returnval.push_back(bkey);
3471 
3472  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussRadauMAlpha1Beta0);
3473  LibUtilities::BasisKey bkey1(LibUtilities::eOrtho_B, nummodes, pkey1);
3474 
3475  returnval.push_back(bkey1);
3476 
3477  const LibUtilities::PointsKey pkey2(nummodes, LibUtilities::eGaussRadauMAlpha2Beta0);
3478  LibUtilities::BasisKey bkey2(LibUtilities::eOrtho_C, nummodes, pkey2);
3479  }
3480  break;
3481  default:
3482  {
3483  ASSERTL0(false,"Expansion not defined in switch for this shape");
3484  }
3485  break;
3486  }
3487  }
3488  break;
3489 
3490  case eGLL_Lagrange_SEM:
3491  {
3492  switch (shape)
3493  {
3494  case LibUtilities::eSegment:
3495  {
3496  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussLobattoLegendre);
3497  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3498 
3499  returnval.push_back(bkey);
3500  }
3501  break;
3502  case LibUtilities::eQuadrilateral:
3503  {
3504  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussLobattoLegendre);
3505  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3506 
3507  returnval.push_back(bkey);
3508  returnval.push_back(bkey);
3509  }
3510  break;
3511  case LibUtilities::eHexahedron:
3512  {
3513  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussLobattoLegendre);
3514  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3515 
3516  returnval.push_back(bkey);
3517  returnval.push_back(bkey);
3518  returnval.push_back(bkey);
3519  }
3520  break;
3521  default:
3522  {
3523  ASSERTL0(false,"Expansion not defined in switch for this shape");
3524  }
3525  break;
3526  }
3527  }
3528  break;
3529 
3530 
3531  case eFourier:
3532  {
3533  switch (shape)
3534  {
3535  case LibUtilities::eSegment:
3536  {
3537  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3538  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3539  returnval.push_back(bkey);
3540  }
3541  break;
3542  case LibUtilities::eQuadrilateral:
3543  {
3544  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3545  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3546  returnval.push_back(bkey);
3547  returnval.push_back(bkey);
3548  }
3549  break;
3550  case LibUtilities::eHexahedron:
3551  {
3552  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3553  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3554  returnval.push_back(bkey);
3555  returnval.push_back(bkey);
3556  returnval.push_back(bkey);
3557  }
3558  break;
3559  default:
3560  {
3561  ASSERTL0(false,"Expansion not defined in switch for this shape");
3562  }
3563  break;
3564  }
3565  }
3566  break;
3567 
3568 
3569  case eFourierSingleMode:
3570  {
3571  switch (shape)
3572  {
3573  case LibUtilities::eSegment:
3574  {
3575  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3576  LibUtilities::BasisKey bkey(LibUtilities::eFourierSingleMode, nummodes, pkey);
3577  returnval.push_back(bkey);
3578  }
3579  break;
3580  case LibUtilities::eQuadrilateral:
3581  {
3582  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3583  LibUtilities::BasisKey bkey(LibUtilities::eFourierSingleMode, nummodes, pkey);
3584  returnval.push_back(bkey);
3585  returnval.push_back(bkey);
3586  }
3587  break;
3588  case LibUtilities::eHexahedron:
3589  {
3590  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3591  LibUtilities::BasisKey bkey(LibUtilities::eFourierSingleMode, nummodes, pkey);
3592  returnval.push_back(bkey);
3593  returnval.push_back(bkey);
3594  returnval.push_back(bkey);
3595  }
3596  break;
3597  default:
3598  {
3599  ASSERTL0(false,"Expansion not defined in switch for this shape");
3600  }
3601  break;
3602  }
3603  }
3604  break;
3605 
3606  case eFourierHalfModeRe:
3607  {
3608  switch (shape)
3609  {
3610  case LibUtilities::eSegment:
3611  {
3612  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3613  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeRe, nummodes, pkey);
3614  returnval.push_back(bkey);
3615  }
3616  break;
3617  case LibUtilities::eQuadrilateral:
3618  {
3619  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3620  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeRe, nummodes, pkey);
3621  returnval.push_back(bkey);
3622  returnval.push_back(bkey);
3623  }
3624  break;
3625  case LibUtilities::eHexahedron:
3626  {
3627  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3628  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeRe, nummodes, pkey);
3629  returnval.push_back(bkey);
3630  returnval.push_back(bkey);
3631  returnval.push_back(bkey);
3632  }
3633  break;
3634  default:
3635  {
3636  ASSERTL0(false,"Expansion not defined in switch for this shape");
3637  }
3638  break;
3639  }
3640  }
3641  break;
3642 
3643  case eFourierHalfModeIm:
3644  {
3645  switch (shape)
3646  {
3647  case LibUtilities::eSegment:
3648  {
3649  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3650  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeIm, nummodes, pkey);
3651  returnval.push_back(bkey);
3652  }
3653  break;
3654  case LibUtilities::eQuadrilateral:
3655  {
3656  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3657  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeIm, nummodes, pkey);
3658  returnval.push_back(bkey);
3659  returnval.push_back(bkey);
3660  }
3661  break;
3662  case LibUtilities::eHexahedron:
3663  {
3664  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3665  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeIm, nummodes, pkey);
3666  returnval.push_back(bkey);
3667  returnval.push_back(bkey);
3668  returnval.push_back(bkey);
3669  }
3670  break;
3671  default:
3672  {
3673  ASSERTL0(false,"Expansion not defined in switch for this shape");
3674  }
3675  break;
3676  }
3677  }
3678  break;
3679 
3680  case eChebyshev:
3681  {
3682  switch (shape)
3683  {
3684  case LibUtilities::eSegment:
3685  {
3686  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussChebyshev);
3687  LibUtilities::BasisKey bkey(LibUtilities::eChebyshev, nummodes, pkey);
3688  returnval.push_back(bkey);
3689  }
3690  break;
3691  case LibUtilities::eQuadrilateral:
3692  {
3693  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussChebyshev);
3694  LibUtilities::BasisKey bkey(LibUtilities::eChebyshev, nummodes, pkey);
3695  returnval.push_back(bkey);
3696  returnval.push_back(bkey);
3697  }
3698  break;
3699  case LibUtilities::eHexahedron:
3700  {
3701  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussChebyshev);
3702  LibUtilities::BasisKey bkey(LibUtilities::eChebyshev, nummodes, pkey);
3703  returnval.push_back(bkey);
3704  returnval.push_back(bkey);
3705  returnval.push_back(bkey);
3706  }
3707  break;
3708  default:
3709  {
3710  ASSERTL0(false,"Expansion not defined in switch for this shape");
3711  }
3712  break;
3713  }
3714  }
3715  break;
3716 
3717  case eFourierChebyshev:
3718  {
3719  switch (shape)
3720  {
3721  case LibUtilities::eQuadrilateral:
3722  {
3723  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3724  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3725  returnval.push_back(bkey);
3726 
3727  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussGaussChebyshev);
3728  LibUtilities::BasisKey bkey1(LibUtilities::eChebyshev, nummodes, pkey1);
3729  returnval.push_back(bkey1);
3730  }
3731  break;
3732  default:
3733  {
3734  ASSERTL0(false,"Expansion not defined in switch for this shape");
3735  }
3736  break;
3737  }
3738  }
3739  break;
3740 
3741  case eChebyshevFourier:
3742  {
3743  switch (shape)
3744  {
3745  case LibUtilities::eQuadrilateral:
3746  {
3747  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussGaussChebyshev);
3748  LibUtilities::BasisKey bkey1(LibUtilities::eChebyshev, nummodes, pkey1);
3749  returnval.push_back(bkey1);
3750 
3751  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3752  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3753  returnval.push_back(bkey);
3754  }
3755  break;
3756  default:
3757  {
3758  ASSERTL0(false,"Expansion not defined in switch for this shape");
3759  }
3760  break;
3761  }
3762  }
3763  break;
3764 
3765  case eFourierModified:
3766  {
3767  switch (shape)
3768  {
3769  case LibUtilities::eQuadrilateral:
3770  {
3771  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3772  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3773  returnval.push_back(bkey);
3774 
3775  const LibUtilities::PointsKey pkey1(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3776  LibUtilities::BasisKey bkey1(LibUtilities::eModified_A, nummodes, pkey1);
3777  returnval.push_back(bkey1);
3778  }
3779  break;
3780  default:
3781  {
3782  ASSERTL0(false,"Expansion not defined in switch for this shape");
3783  }
3784  break;
3785  }
3786  }
3787  break;
3788 
3789  default:
3790  {
3791  ASSERTL0(false,"Expansion type not defined");
3792  }
3793  break;
3794 
3795  }
3796 
3797  return returnval;
3798  }
3799 
3800  /**
3801  *
3802  */
3803  LibUtilities::BasisKeyVector
3804  MeshGraph::DefineBasisKeyFromExpansionTypeHomo(
3805  GeometrySharedPtr in,
3806  ExpansionType type_x,
3807  ExpansionType type_y,
3808  ExpansionType type_z,
3809  const int nummodes_x,
3810  const int nummodes_y,
3811  const int nummodes_z)
3812  {
3813  LibUtilities::BasisKeyVector returnval;
3814 
3815  LibUtilities::ShapeType shape = in->GetShapeType();
3816 
3817  switch (shape)
3818  {
3819  case LibUtilities::eSegment:
3820  {
3821  ASSERTL0(false,"Homogeneous expansion not defined for this shape");
3822  }
3823  break;
3824 
3825  case LibUtilities::eQuadrilateral:
3826  {
3827  ASSERTL0(false,"Homogeneous expansion not defined for this shape");
3828  }
3829  break;
3830 
3831  case LibUtilities::eHexahedron:
3832  {
3833  switch(type_x)
3834  {
3835  case eFourier:
3836  {
3837  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eFourierEvenlySpaced);
3838  LibUtilities::BasisKey bkey1(LibUtilities::eFourier,nummodes_x,pkey1);
3839  returnval.push_back(bkey1);
3840  }
3841  break;
3842 
3843  case eFourierSingleMode:
3844  {
3845  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eFourierSingleModeSpaced);
3846  LibUtilities::BasisKey bkey1(LibUtilities::eFourierSingleMode,nummodes_x,pkey1);
3847  returnval.push_back(bkey1);
3848  }
3849  break;
3850 
3851  case eFourierHalfModeRe:
3852  {
3853  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eFourierSingleModeSpaced);
3854  LibUtilities::BasisKey bkey1(LibUtilities::eFourierHalfModeRe,nummodes_x,pkey1);
3855  returnval.push_back(bkey1);
3856  }
3857  break;
3858 
3859  case eFourierHalfModeIm:
3860  {
3861  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eFourierSingleModeSpaced);
3862  LibUtilities::BasisKey bkey1(LibUtilities::eFourierHalfModeIm,nummodes_x,pkey1);
3863  returnval.push_back(bkey1);
3864  }
3865  break;
3866 
3867 
3868  case eChebyshev:
3869  {
3870  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eGaussGaussChebyshev);
3871  LibUtilities::BasisKey bkey1(LibUtilities::eChebyshev,nummodes_x,pkey1);
3872  returnval.push_back(bkey1);
3873  }
3874  break;
3875 
3876 
3877 
3878  default:
3879  {
3880  ASSERTL0(false,"Homogeneous expansion can be of Fourier or Chebyshev type only");
3881  }
3882  break;
3883  }
3884 
3885 
3886  switch(type_y)
3887  {
3888  case eFourier:
3889  {
3890  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eFourierEvenlySpaced);
3891  LibUtilities::BasisKey bkey2(LibUtilities::eFourier,nummodes_y,pkey2);
3892  returnval.push_back(bkey2);
3893  }
3894  break;
3895 
3896 
3897  case eFourierSingleMode:
3898  {
3899  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eFourierSingleModeSpaced);
3900  LibUtilities::BasisKey bkey2(LibUtilities::eFourierSingleMode,nummodes_y,pkey2);
3901  returnval.push_back(bkey2);
3902  }
3903  break;
3904 
3905  case eFourierHalfModeRe:
3906  {
3907  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eFourierSingleModeSpaced);
3908  LibUtilities::BasisKey bkey2(LibUtilities::eFourierHalfModeRe,nummodes_y,pkey2);
3909  returnval.push_back(bkey2);
3910  }
3911  break;
3912 
3913  case eFourierHalfModeIm:
3914  {
3915  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eFourierSingleModeSpaced);
3916  LibUtilities::BasisKey bkey2(LibUtilities::eFourierHalfModeIm,nummodes_y,pkey2);
3917  returnval.push_back(bkey2);
3918  }
3919  break;
3920 
3921  case eChebyshev:
3922  {
3923  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eGaussGaussChebyshev);
3924  LibUtilities::BasisKey bkey2(LibUtilities::eChebyshev,nummodes_y,pkey2);
3925  returnval.push_back(bkey2);
3926  }
3927  break;
3928 
3929  default:
3930  {
3931  ASSERTL0(false,"Homogeneous expansion can be of Fourier or Chebyshev type only");
3932  }
3933  break;
3934  }
3935 
3936  switch(type_z)
3937  {
3938  case eFourier:
3939  {
3940  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eFourierEvenlySpaced);
3941  LibUtilities::BasisKey bkey3(LibUtilities::eFourier,nummodes_z,pkey3);
3942  returnval.push_back(bkey3);
3943  }
3944  break;
3945 
3946  case eFourierSingleMode:
3947  {
3948  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eFourierSingleModeSpaced);
3949  LibUtilities::BasisKey bkey3(LibUtilities::eFourierSingleMode,nummodes_z,pkey3);
3950  returnval.push_back(bkey3);
3951  }
3952  break;
3953 
3954  case eFourierHalfModeRe:
3955  {
3956  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eFourierSingleModeSpaced);
3957  LibUtilities::BasisKey bkey3(LibUtilities::eFourierHalfModeRe,nummodes_z,pkey3);
3958  returnval.push_back(bkey3);
3959  }
3960  break;
3961 
3962  case eFourierHalfModeIm:
3963  {
3964  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eFourierSingleModeSpaced);
3965  LibUtilities::BasisKey bkey3(LibUtilities::eFourierHalfModeIm,nummodes_z,pkey3);
3966  returnval.push_back(bkey3);
3967  }
3968  break;
3969 
3970  case eChebyshev:
3971  {
3972  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eGaussGaussChebyshev);
3973  LibUtilities::BasisKey bkey3(LibUtilities::eChebyshev,nummodes_z,pkey3);
3974  returnval.push_back(bkey3);
3975  }
3976  break;
3977 
3978  default:
3979  {
3980  ASSERTL0(false,"Homogeneous expansion can be of Fourier or Chebyshev type only");
3981  }
3982  break;
3983  }
3984  }
3985  break;
3986 
3987  case LibUtilities::eTriangle:
3988  {
3989  ASSERTL0(false,"Homogeneous expansion not defined for this shape");
3990  }
3991  break;
3992 
3993  case LibUtilities::eTetrahedron:
3994  {
3995  ASSERTL0(false,"Homogeneous expansion not defined for this shape");
3996  }
3997  break;
3998 
3999  default:
4000  ASSERTL0(false,"Expansion not defined in switch for this shape");
4001  break;
4002  }
4003 
4004  return returnval;
4005  }
4006 
4007 
4008  /**
4009  *
4010  */
4011  PointGeomSharedPtr MeshGraph::AddVertex(NekDouble x, NekDouble y, NekDouble z)
4012  {
4013  unsigned int nextId = m_vertSet.rbegin()->first + 1;
4014  PointGeomSharedPtr vert(MemoryManager<PointGeom>::AllocateSharedPtr(m_spaceDimension, nextId, x, y, z));
4015  m_vertSet[nextId] = vert;
4016  return vert;
4017  }
4018 
4019  /**
4020  *
4021  */
4022  SegGeomSharedPtr MeshGraph::AddEdge(PointGeomSharedPtr v0, PointGeomSharedPtr v1,
4023  CurveSharedPtr curveDefinition)
4024  {
4025  PointGeomSharedPtr vertices[] = {v0, v1};
4026  SegGeomSharedPtr edge;
4027  int edgeId = m_segGeoms.rbegin()->first + 1;
4028 
4029  if( curveDefinition )
4030  {
4031  edge = MemoryManager<SegGeom>::AllocateSharedPtr(edgeId, m_spaceDimension, vertices, curveDefinition);
4032  }
4033  else
4034  {
4035  edge = MemoryManager<SegGeom>::AllocateSharedPtr(edgeId, m_spaceDimension, vertices);
4036  }
4037  m_segGeoms[edgeId] = edge;
4038  return edge;
4039  }
4040 
4041 
4042  /**
4043  *
4044  */
4045  TriGeomSharedPtr MeshGraph::AddTriangle(SegGeomSharedPtr edges[], StdRegions::Orientation orient[])
4046  {
4047  int indx = m_triGeoms.rbegin()->first + 1;
4048  TriGeomSharedPtr trigeom(MemoryManager<TriGeom>::AllocateSharedPtr(indx, edges, orient));
4049  trigeom->SetGlobalID(indx);
4050 
4051  m_triGeoms[indx] = trigeom;
4052 
4053  return trigeom;
4054  }
4055 
4056 
4057  /**
4058  *
4059  */
4060  QuadGeomSharedPtr MeshGraph::AddQuadrilateral(SegGeomSharedPtr edges[], StdRegions::Orientation orient[])
4061  {
4062  int indx = m_quadGeoms.rbegin()->first + 1;
4063  QuadGeomSharedPtr quadgeom(MemoryManager<QuadGeom>::AllocateSharedPtr(indx, edges, orient));
4064  quadgeom->SetGlobalID(indx);
4065 
4066  m_quadGeoms[indx] = quadgeom;
4067  return quadgeom;
4068  }
4069 
4070 
4071  /**
4072  *
4073  */
4074  PrismGeomSharedPtr MeshGraph::AddPrism(TriGeomSharedPtr tfaces[PrismGeom::kNtfaces],
4075  QuadGeomSharedPtr qfaces[PrismGeom::kNqfaces])
4076  {
4077  // Setting the orientation is disabled in the reader. Why?
4078  Geometry2DSharedPtr faces[] = { qfaces[0], tfaces[0], qfaces[1], tfaces[1], qfaces[2] };
4079  unsigned int index = m_prismGeoms.rbegin()->first + 1;
4080  PrismGeomSharedPtr prismgeom(MemoryManager<PrismGeom>::AllocateSharedPtr(faces));
4081  prismgeom->SetGlobalID(index);
4082 
4083  m_prismGeoms[index] = prismgeom;
4084  return prismgeom;
4085  }
4086 
4087 
4088  /**
4089  *
4090  */
4091  TetGeomSharedPtr MeshGraph::AddTetrahedron(TriGeomSharedPtr tfaces[TetGeom::kNtfaces])
4092  {
4093  unsigned int index = m_tetGeoms.rbegin()->first + 1;
4094  TetGeomSharedPtr tetgeom(MemoryManager<TetGeom>::AllocateSharedPtr(tfaces));
4095  tetgeom->SetGlobalID(index);
4096 
4097  m_tetGeoms[index] = tetgeom;
4098  return tetgeom;
4099  }
4100 
4101 
4102  /**
4103  *
4104  */
4105  PyrGeomSharedPtr MeshGraph::AddPyramid(TriGeomSharedPtr tfaces[PyrGeom::kNtfaces],
4106  QuadGeomSharedPtr qfaces[PyrGeom::kNqfaces])
4107  {
4108  Geometry2DSharedPtr faces[] = { qfaces[0], tfaces[0], tfaces[1], tfaces[2], tfaces[3] };
4109  unsigned int index = m_pyrGeoms.rbegin()->first + 1;
4110 
4111  PyrGeomSharedPtr pyrgeom(MemoryManager<PyrGeom>::AllocateSharedPtr(faces));
4112  pyrgeom->SetGlobalID(index);
4113 
4114  m_pyrGeoms[index] = pyrgeom;
4115  return pyrgeom;
4116  }
4117 
4118 
4119  /**
4120  *
4121  */
4122  HexGeomSharedPtr MeshGraph::AddHexahedron(QuadGeomSharedPtr qfaces[HexGeom::kNqfaces])
4123  {
4124  unsigned int index = m_hexGeoms.rbegin()->first + 1;
4125  HexGeomSharedPtr hexgeom(MemoryManager<HexGeom>::AllocateSharedPtr(qfaces));
4126  hexgeom->SetGlobalID(index);
4127  m_hexGeoms[index] = hexgeom;
4128  return hexgeom;
4129  }
4130 
4131 
4132  /**
4133  * Generate a single vector of Expansion structs mapping global element
4134  * ID to a corresponding Geometry shared pointer and basis key.
4135  *
4136  * Expansion map ensures elements which appear in multiple composites
4137  * within the domain are only listed once.
4138  */
4139  ExpansionMapShPtr MeshGraph::SetUpExpansionMap(void)
4140  {
4141  ExpansionMapShPtr returnval;
4142  returnval = MemoryManager<ExpansionMap>::AllocateSharedPtr();
4143 
4144  for(int d = 0; d < m_domain.size(); ++d)
4145  {
4146  CompositeMap::const_iterator compIter;
4147 
4148  for (compIter = m_domain[d].begin(); compIter != m_domain[d].end(); ++compIter)
4149  {
4150  GeometryVector::const_iterator x;
4151  for (x = compIter->second->begin(); x != compIter->second->end(); ++x)
4152  {
4153  LibUtilities::BasisKeyVector def;
4154  ExpansionShPtr expansionElementShPtr =
4155  MemoryManager<Expansion>::AllocateSharedPtr(*x, def);
4156  int id = (*x)->GetGlobalID();
4157  (*returnval)[id] = expansionElementShPtr;
4158  }
4159  }
4160  }
4161 
4162  return returnval;
4163  }
4164  }; //end of namespace
4165 }; //end of namespace
StdPrismExp.h
Nektar::SpatialDomains::MeshGraph::SetExpansions
void SetExpansions(std::vector< LibUtilities::FieldDefinitionsSharedPtr > &fielddef)
Sets expansions given field definitions.
Definition: MeshGraph.cpp:2369
Nektar::SpatialDomains::PyrGeomSharedPtr
boost::shared_ptr< PyrGeom > PyrGeomSharedPtr
Definition: PyrGeom.h:84
Nektar::NekInt
boost::int32_t NekInt
Definition: NektarUnivTypeDefs.hpp:46
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Nektar::ParseUtils::GenerateOrderedStringVector
static bool GenerateOrderedStringVector(const char *const str, std::vector< std::string > &vec)
Definition: ParseUtils.hpp:143
Nektar::SpatialDomains::MeshGraph::m_session
LibUtilities::SessionReaderSharedPtr m_session
Definition: MeshGraph.h:409
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:191
Nektar::LibUtilities::eModified_C
Principle Modified Functions .
Definition: BasisType.h:51
Nektar::ParseUtils::GenerateOrderedVector
static bool GenerateOrderedVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:97
Nektar::SpatialDomains::MeshGraph::SetExpansionsToPolyOrder
void SetExpansionsToPolyOrder(int nmodes)
Reset expansion to have specified polynomial order nmodes.
Definition: MeshGraph.cpp:3111
Nektar::SpatialDomains::MeshGraph::Read
static boost::shared_ptr< MeshGraph > Read(const LibUtilities::SessionReaderSharedPtr &pSession, DomainRangeShPtr &rng=NullDomainRangeShPtr)
Definition: MeshGraph.cpp:124
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Nektar::LibUtilities::BasisTypeMap
const char *const BasisTypeMap[]
Definition: Foundations.hpp:47
Nektar::SpatialDomains::MeshGraph::ReadCurves
void ReadCurves(TiXmlDocument &doc)
Definition: MeshGraph.cpp:1181
Nektar::SpatialDomains::MeshGraph::AddPrism
PrismGeomSharedPtr AddPrism(TriGeomSharedPtr tfaces[PrismGeom::kNtfaces], QuadGeomSharedPtr qfaces[PrismGeom::kNqfaces])
Definition: MeshGraph.cpp:4074
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:102
Nektar::LibUtilities::kPointsTypeStr
const std::string kPointsTypeStr[]
Definition: Foundations.hpp:69
FieldIOXml.h
Nektar::LibUtilities::BasisKey::GetBasisType
BasisType GetBasisType() const
Return type of expansion basis.
Definition: Basis.h:139
Nektar::SpatialDomains::QuadGeomSharedPtr
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: QuadGeom.h:54
Nektar::LibUtilities::MeshCurvedPts::index
std::vector< NekInt64 > index
id of this Point set
Definition: MeshEntities.hpp:130
Nektar::SpatialDomains::MeshGraph::GetMeshDimension
int GetMeshDimension() const
Dimension of the mesh (can be a 1D curve in 3D space).
Definition: MeshGraph.h:448
Nektar::LibUtilities::eModified_A
Principle Modified Functions .
Definition: BasisType.h:49
Nektar::LibUtilities::eGauss_Lagrange
Lagrange Polynomials using the Gauss points .
Definition: BasisType.h:54
Nektar::ParseUtils::GenerateSeqString
static std::string GenerateSeqString(const std::vector< unsigned int > &elmtids)
Definition: ParseUtils.hpp:159
Nektar::LibUtilities::AnalyticExpressionEvaluator::Evaluate
NekDouble Evaluate(const int AnalyticExpression_id)
Evaluation method for expressions depending on parameters only.
Definition: AnalyticExpressionEvaluator.cpp:494
std
STL namespace.
Nektar::LibUtilities::BasisKey::GetPointsType
PointsType GetPointsType() const
Return type of quadrature.
Definition: Basis.h:151
Nektar::LibUtilities::eFourier
Fourier Expansion .
Definition: BasisType.h:52
Nektar::LibUtilities::eChebyshev
Chebyshev Polynomials .
Definition: BasisType.h:56
Nektar::SpatialDomains::TetGeom::kNtfaces
static const int kNtfaces
Definition: TetGeom.h:59
Nektar::SpatialDomains::HexGeomSharedPtr
boost::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:110
Nektar::SpatialDomains::MeshGraph::AddVertex
PointGeomSharedPtr AddVertex(NekDouble x, NekDouble y, NekDouble z)
Adds a vertex to the with the next available ID.
Definition: MeshGraph.cpp:4011
Nektar::SpatialDomains::MeshGraph::AddQuadrilateral
QuadGeomSharedPtr AddQuadrilateral(SegGeomSharedPtr edges[], StdRegions::Orientation orient[])
Definition: MeshGraph.cpp:4060
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Nektar::ParseUtils::GenerateSeqVector
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:79
Nektar::SpatialDomains::MeshGraph::ReadGeometry
virtual void ReadGeometry(const std::string &infilename)
Read will read the meshgraph vertices given a filename.
Definition: MeshGraph.cpp:235
MeshGraph1D.h
Nektar::SpatialDomains::CurveSharedPtr
boost::shared_ptr< Curve > CurveSharedPtr
Definition: Curve.hpp:62
Nektar::SpatialDomains::MeshGraph::GetExpansions
const ExpansionMap & GetExpansions()
Definition: MeshGraph.h:515
Nektar::LibUtilities::eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:49
Nektar::LibUtilities::eGaussRadauMAlpha1Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:59
Nektar::SpatialDomains::MeshGraph::DefineBasisKeyFromExpansionType
static LibUtilities::BasisKeyVector DefineBasisKeyFromExpansionType(GeometrySharedPtr in, ExpansionType type, const int order)
Definition: MeshGraph.cpp:3212
CellMLToNektar.cellml_metadata.p
tuple p
Definition: cellml_metadata.py:350
Nektar::SpatialDomains::GeometryVectorIter
std::vector< GeometrySharedPtr >::iterator GeometryVectorIter
Definition: Geometry.h:58
Nektar::LibUtilities::eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:47
Nektar::SpatialDomains::MeshGraph::SetExpansionsToPointOrder
void SetExpansionsToPointOrder(int npts)
Reset expansion to have specified point order npts.
Definition: MeshGraph.cpp:3144
Nektar::LibUtilities::ePolyEvenlySpaced
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:65
Nektar::SpatialDomains::MeshGraph::AddTriangle
TriGeomSharedPtr AddTriangle(SegGeomSharedPtr edges[], StdRegions::Orientation orient[])
Definition: MeshGraph.cpp:4045
Nektar::SpatialDomains::MeshGraph::m_domainRange
DomainRangeShPtr m_domainRange
Definition: MeshGraph.h:433
Nektar::LibUtilities::eFourierEvenlySpaced
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:66
Nektar::LibUtilities::eFourierHalfModeRe
Fourier Modified expansions with just the real part of the first mode .
Definition: BasisType.h:59
Nektar::SpatialDomains::MeshGraph::WriteGeometry
void WriteGeometry(std::string &outfilename)
Write out an XML file containing the GEOMETRY block representing this MeshGraph instance inside a NEK...
Definition: MeshGraph.cpp:1977
Nektar::SpatialDomains::MeshGraph::GetCompositeItem
GeometrySharedPtr GetCompositeItem(int whichComposite, int whichItem)
Definition: MeshGraph.cpp:2239
Nektar::SpatialDomains::ExpansionMapShPtr
boost::shared_ptr< ExpansionMap > ExpansionMapShPtr
Definition: MeshGraph.h:178
Nektar::SpatialDomains::MeshGraph::AddEdge
SegGeomSharedPtr AddEdge(PointGeomSharedPtr v0, PointGeomSharedPtr v1, CurveSharedPtr curveDefinition=CurveSharedPtr())
Adds an edge between two points. If curveDefinition is null, then the edge is straight, otherwise it is curved according to the curveDefinition.
Definition: MeshGraph.cpp:4022
Nektar::SpatialDomains::MeshGraph::GetComposite
Composite GetComposite(int whichComposite) const
Definition: MeshGraph.h:466
Nektar::SpatialDomains::SegGeomSharedPtr
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
Nektar::LibUtilities::eModified_B
Principle Modified Functions .
Definition: BasisType.h:50
Nektar::LibUtilities::XmlDataSource::create
static DataSourceSharedPtr create(const std::string &fn)
Create a new XML data source based on the filename.
Definition: FieldIOXml.h:91
Nektar::SpatialDomains::ExpansionMapShPtrMapIter
std::map< std::string, ExpansionMapShPtr >::iterator ExpansionMapShPtrMapIter
Definition: MeshGraph.h:180
StdPyrExp.h
Nektar::LibUtilities::Equation::Evaluate
NekDouble Evaluate() const
Definition: Equation.h:102
Nektar::LibUtilities::BasisKey::GetNumPoints
int GetNumPoints() const
Return points order at which basis is defined.
Definition: Basis.h:128
npts
static std::string npts
Definition: InputFld.cpp:43
Nektar::SpatialDomains::DomainRangeShPtr
boost::shared_ptr< DomainRange > DomainRangeShPtr
Definition: MeshGraph.h:157
Nektar::SpatialDomains::MeshGraph::AddPyramid
PyrGeomSharedPtr AddPyramid(TriGeomSharedPtr tfaces[PyrGeom::kNtfaces], QuadGeomSharedPtr qfaces[PyrGeom::kNqfaces])
Definition: MeshGraph.cpp:4105
Nektar::LibUtilities::eOrtho_C
Principle Orthogonal Functions .
Definition: BasisType.h:48
Nektar::SpatialDomains::CompositeMapIter
std::map< int, Composite >::iterator CompositeMapIter
Definition: MeshGraph.h:116
Nektar::LibUtilities::eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:46
Nektar::LibUtilities::eGaussGaussChebyshev
1D Gauss-Gauss-Chebyshev quadrature points
Definition: PointsType.h:53
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:58
Nektar::LibUtilities::FieldIOSharedPtr
boost::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:309
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::SpatialDomains::MeshGraph::DefineBasisKeyFromExpansionTypeHomo
LibUtilities::BasisKeyVector DefineBasisKeyFromExpansionTypeHomo(GeometrySharedPtr in, ExpansionType type_x, ExpansionType type_y, ExpansionType type_z, const int nummodes_x, const int nummodes_y, const int nummodes_z)
Definition: MeshGraph.cpp:3804
StdTriExp.h
Nektar::LibUtilities::BasisKeyVector
std::vector< BasisKey > BasisKeyVector
Name for a vector of BasisKeys.
Definition: FoundationsFwd.hpp:70
MeshGraph3D.h
Nektar::SpatialDomains::Composite
boost::shared_ptr< GeometryVector > Composite
Definition: MeshGraph.h:114
Nektar::NekConstants::kNekUnsetDouble
static const NekDouble kNekUnsetDouble
Definition: NektarUnivConsts.hpp:47
Nektar::LibUtilities::AnalyticExpressionEvaluator
This class defines evaluator of analytic (symbolic) mathematical expressions. Expressions are allowed...
Definition: AnalyticExpressionEvaluator.hpp:114
Nektar::SpatialDomains::Geometry::GetVertex
PointGeomSharedPtr GetVertex(int i) const
Definition: Geometry.h:348
Nektar::SpatialDomains::CompositeMap
std::map< int, Composite > CompositeMap
Definition: MeshGraph.h:115
Nektar::LibUtilities::eFourierHalfModeIm
Fourier Modified expansions with just the imaginary part of the first mode .
Definition: BasisType.h:60
Nektar::SpatialDomains::PyrGeom::kNtfaces
static const int kNtfaces
Definition: PyrGeom.h:59
MeshGraph2D.h
Nektar::SpatialDomains::Geometry2DSharedPtr
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry2D.h:59
Nektar::SpatialDomains::MeshGraph::ReadGeometryInfo
void ReadGeometryInfo(const std::string &infilename)
Read geometric information from a file.
Definition: MeshGraph.cpp:530
Nektar::iterator
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
Definition: MatrixOperations.cpp:96
MeshGraph.h
Nektar::LibUtilities::AnalyticExpressionEvaluator::DefineFunction
int DefineFunction(const std::string &vlist, const std::string &function)
This function allows one to define a function to evaluate. The first argument (vlist) is a list of va...
Definition: AnalyticExpressionEvaluator.cpp:402
Nektar::LibUtilities::FieldIO::CreateForFile
static boost::shared_ptr< FieldIO > CreateForFile(const LibUtilities::SessionReaderSharedPtr session, const std::string &filename)
Construct a FieldIO object for a given input filename.
Definition: FieldIO.cpp:212
Nektar::SpatialDomains::PyrGeom::kNqfaces
static const int kNqfaces
Definition: PyrGeom.h:58
Nektar::SpatialDomains::Geometry3D
3D geometry information
Definition: Geometry3D.h:70
Nektar::SpatialDomains::MeshGraph::m_domain
std::vector< CompositeMap > m_domain
Definition: MeshGraph.h:432
Nektar::SpatialDomains::MeshGraph::GetCompositeList
void GetCompositeList(const std::string &compositeStr, CompositeMap &compositeVector) const
Definition: MeshGraph.cpp:2277
Nektar::SpatialDomains::NullDomainRangeShPtr
static DomainRangeShPtr NullDomainRangeShPtr
Definition: MeshGraph.h:158
Nektar::SpatialDomains::Geometry2D
2D geometry information
Definition: Geometry2D.h:65
Nektar::SpatialDomains::ExpansionShPtr
boost::shared_ptr< Expansion > ExpansionShPtr
Definition: MeshGraph.h:173
Nektar::SpatialDomains::MeshGraph::m_expansionMapShPtrMap
ExpansionMapShPtrMap m_expansionMapShPtrMap
Definition: MeshGraph.h:435
Nektar::LibUtilities::eFourierSingleMode
Fourier ModifiedExpansion with just the first mode .
Definition: BasisType.h:58
Nektar::SpatialDomains::MeshGraph
Base class for a spectral/hp element mesh.
Definition: MeshGraph.h:186
Equation.h
Nektar::SpatialDomains::MeshGraph::ReadDomain
void ReadDomain(TiXmlDocument &doc)
Definition: MeshGraph.cpp:1087
Nektar::SpatialDomains::PrismGeomSharedPtr
boost::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:109
Nektar::SpatialDomains::MeshGraph::SetDomainRange
void SetDomainRange(NekDouble xmin, NekDouble xmax, NekDouble ymin=NekConstants::kNekUnsetDouble, NekDouble ymax=NekConstants::kNekUnsetDouble, NekDouble zmin=NekConstants::kNekUnsetDouble, NekDouble zmax=NekConstants::kNekUnsetDouble)
Definition: MeshGraph.cpp:1991
Nektar::LibUtilities::eFourierSingleModeSpaced
1D Non Evenly-spaced points for Single Mode analysis
Definition: PointsType.h:67
Nektar::SpatialDomains::MeshGraph::CheckRange
bool CheckRange(Geometry2D &geom)
Check if goemetry is in range definition if activated.
Definition: MeshGraph.cpp:2029
Nektar::LibUtilities::eGaussRadauMAlpha2Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:60
Nektar::LibUtilities::MeshCurvedPts::pts
std::vector< MeshVertex > pts
mapping to access pts value.
Definition: MeshEntities.hpp:134
Nektar::SpatialDomains::MeshGraph::SetBasisKey
void SetBasisKey(LibUtilities::ShapeType shape, LibUtilities::BasisKeyVector &keys, std::string var="DefaultVar")
Sets the basis key for all expansions of the given shape.
Definition: MeshGraph.cpp:3191
Nektar::SpatialDomains::MeshGraph::GetExpansion
ExpansionShPtr GetExpansion(GeometrySharedPtr geom, const std::string variable="DefaultVar")
Definition: MeshGraph.cpp:2352
Nektar::SpatialDomains::TetGeomSharedPtr
boost::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:106
Nektar::SpatialDomains::TriGeomSharedPtr
boost::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
Nektar::SpatialDomains::MeshGraph::AddHexahedron
HexGeomSharedPtr AddHexahedron(QuadGeomSharedPtr qfaces[HexGeom::kNqfaces])
Definition: MeshGraph.cpp:4122
Nektar::StdRegions::find
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:316
Nektar::SpatialDomains::MeshGraph::SetUpExpansionMap
ExpansionMapShPtr SetUpExpansionMap(void)
Definition: MeshGraph.cpp:4139
Nektar::SpatialDomains::MeshGraph::ReadExpansions
void ReadExpansions(const std::string &infilename)
Read the expansions given the XML file path.
Definition: MeshGraph.cpp:585
Nektar::SpatialDomains::HexGeom::kNqfaces
static const int kNqfaces
Definition: HexGeom.h:62
Nektar::LibUtilities::BasisKey::GetNumModes
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:84
Nektar::SpatialDomains::kExpansionTypeStr
const std::string kExpansionTypeStr[]
Definition: MeshGraph.h:86
Nektar::LibUtilities::eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:53
ParseUtils.hpp
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
CompressData.h
Nektar::SpatialDomains::GeometrySharedPtr
boost::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:53
Nektar::SpatialDomains::MeshGraph::SetExpansionsToEvenlySpacedPoints
void SetExpansionsToEvenlySpacedPoints(int npoints=0)
Sets expansions to have equispaced points.
Definition: MeshGraph.cpp:3069
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
Nektar::SpatialDomains::Geometry::GetShapeType
LibUtilities::ShapeType GetShapeType(void)
Definition: Geometry.h:304
Nektar::SpatialDomains::PointGeomSharedPtr
boost::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60
Nektar::LibUtilities::eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:52
Nektar::SpatialDomains::ExpansionMap
std::map< int, ExpansionShPtr > ExpansionMap
Definition: MeshGraph.h:174
Nektar::LibUtilities::CompressData::ZlibDecodeFromBase64Str
int ZlibDecodeFromBase64Str(std::string &in64, std::vector< T > &out)
Definition: CompressData.h:243
Nektar::SpatialDomains::MeshGraph::AddTetrahedron
TetGeomSharedPtr AddTetrahedron(TriGeomSharedPtr tfaces[TetGeom::kNtfaces])
Definition: MeshGraph.cpp:4091
StdTetExp.h
Nektar::SpatialDomains::ExpansionMapIter
std::map< int, ExpansionShPtr >::iterator ExpansionMapIter
Definition: MeshGraph.h:175
Generated on Sat Nov 25 2017 16:18:14 for Nektar++ by   doxygen 1.8.8