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  GeometrySharedPtr firstGeom = cIt->second->at(0);
1929  int shapeDim = firstGeom->GetShapeDim();
1930  string tag = (shapeDim < m_meshDimension) ?
1931  compMap[firstGeom->GetShapeType()].second :
1932  compMap[firstGeom->GetShapeType()].first;
1933 
1934  idxList.clear();
1935  s << " " << tag << "[";
1936 
1937  for (int i = 0; i < cIt->second->size(); ++i)
1938  {
1939  idxList.push_back((*cIt->second)[i]->GetGlobalID());
1940  }
1941 
1942  s << ParseUtils::GenerateSeqString(idxList) << "] ";
1943 
1944  c->SetAttribute("ID", cIt->first);
1945  c->LinkEndChild(new TiXmlText(s.str()));
1946  compTag->LinkEndChild(c);
1947  }
1948 
1949  geomTag->LinkEndChild(compTag);
1950 
1951  // Construct <DOMAIN> block
1952  TiXmlElement *domTag = new TiXmlElement("DOMAIN");
1953  stringstream domString;
1954 
1955  // @todo Fix this to accomodate multi domain output
1956  idxList.clear();
1957  for (cIt = m_domain[0].begin(); cIt != m_domain[0].end(); ++cIt)
1958  {
1959  idxList.push_back(cIt->first);
1960  }
1961 
1962  domString << " C[" << ParseUtils::GenerateSeqString(idxList) << "] ";
1963  domTag->LinkEndChild(new TiXmlText(domString.str()));
1964  geomTag->LinkEndChild(domTag);
1965  }
1966 
1967  /**
1968  * @brief Write out an XML file containing the GEOMETRY block
1969  * representing this MeshGraph instance inside a NEKTAR tag.
1970  */
1971  void MeshGraph::WriteGeometry(std::string &outfilename)
1972  {
1973  // Create empty TinyXML document.
1974  TiXmlDocument doc;
1975  TiXmlDeclaration* decl = new TiXmlDeclaration( "1.0", "utf-8", "");
1976  doc.LinkEndChild(decl);
1977 
1978  // Write out geometry information.
1979  WriteGeometry(doc);
1980 
1981  // Save file.
1982  doc.SaveFile(outfilename);
1983  }
1984 
1985  void MeshGraph::SetDomainRange (
1986  NekDouble xmin, NekDouble xmax, NekDouble ymin,
1987  NekDouble ymax, NekDouble zmin, NekDouble zmax)
1988  {
1989  m_domainRange->m_checkShape = false;
1990 
1991  if(m_domainRange == NullDomainRangeShPtr)
1992  {
1993  m_domainRange = MemoryManager<DomainRange>::AllocateSharedPtr();
1994  m_domainRange->m_doXrange = true;
1995  }
1996 
1997  m_domainRange->m_xmin = xmin;
1998  m_domainRange->m_xmax = xmax;
1999 
2000  if(ymin == NekConstants::kNekUnsetDouble)
2001  {
2002  m_domainRange->m_doYrange = false;
2003  }
2004  else
2005  {
2006  m_domainRange->m_doYrange = true;
2007  m_domainRange->m_ymin = ymin;
2008  m_domainRange->m_ymax = ymax;
2009  }
2010 
2011  if(zmin == NekConstants::kNekUnsetDouble)
2012  {
2013  m_domainRange->m_doZrange = false;
2014  }
2015  else
2016  {
2017  m_domainRange->m_doZrange = true;
2018  m_domainRange->m_zmin = zmin;
2019  m_domainRange->m_zmax = zmax;
2020  }
2021  }
2022 
2023  bool MeshGraph::CheckRange(Geometry2D &geom)
2024  {
2025  bool returnval = true;
2026 
2027  if(m_domainRange != NullDomainRangeShPtr)
2028  {
2029  int nverts = geom.GetNumVerts();
2030  int coordim = geom.GetCoordim();
2031 
2032  // exclude elements outside x range if all vertices not in region
2033  if(m_domainRange->m_doXrange)
2034  {
2035  int ncnt_low = 0;
2036  int ncnt_up = 0;
2037  for(int i = 0; i < nverts; ++i)
2038  {
2039  NekDouble xval = (*geom.GetVertex(i))[0];
2040  if(xval < m_domainRange->m_xmin)
2041  {
2042  ncnt_low++;
2043  }
2044 
2045  if(xval > m_domainRange->m_xmax)
2046  {
2047  ncnt_up++;
2048  }
2049  }
2050 
2051  // check for all verts to be less or greater than
2052  // range so that if element spans thin range then
2053  // it is still included
2054  if((ncnt_up == nverts)||(ncnt_low == nverts))
2055  {
2056  returnval = false;
2057  }
2058  }
2059 
2060  // exclude elements outside y range if all vertices not in region
2061  if(m_domainRange->m_doYrange)
2062  {
2063  int ncnt_low = 0;
2064  int ncnt_up = 0;
2065  for(int i = 0; i < nverts; ++i)
2066  {
2067  NekDouble yval = (*geom.GetVertex(i))[1];
2068  if(yval < m_domainRange->m_ymin)
2069  {
2070  ncnt_low++;
2071  }
2072 
2073  if(yval > m_domainRange->m_ymax)
2074  {
2075  ncnt_up++;
2076  }
2077  }
2078 
2079  // check for all verts to be less or greater than
2080  // range so that if element spans thin range then
2081  // it is still included
2082  if((ncnt_up == nverts)||(ncnt_low == nverts))
2083  {
2084  returnval = false;
2085  }
2086  }
2087 
2088  if(coordim > 2)
2089  {
2090  // exclude elements outside z range if all vertices not in region
2091  if(m_domainRange->m_doZrange)
2092  {
2093  int ncnt_low = 0;
2094  int ncnt_up = 0;
2095 
2096  for(int i = 0; i < nverts; ++i)
2097  {
2098  NekDouble zval = (*geom.GetVertex(i))[2];
2099 
2100  if(zval < m_domainRange->m_zmin)
2101  {
2102  ncnt_low++;
2103  }
2104 
2105  if(zval > m_domainRange->m_zmax)
2106  {
2107  ncnt_up++;
2108  }
2109  }
2110 
2111  // check for all verts to be less or greater than
2112  // range so that if element spans thin range then
2113  // it is still included
2114  if((ncnt_up == nverts)||(ncnt_low == nverts))
2115  {
2116  returnval = false;
2117  }
2118  }
2119  }
2120  }
2121  return returnval;
2122  }
2123 
2124 
2125  /* Domain checker for 3D geometries */
2126  bool MeshGraph::CheckRange(Geometry3D &geom)
2127  {
2128  bool returnval = true;
2129 
2130  if(m_domainRange != NullDomainRangeShPtr)
2131  {
2132  int nverts = geom.GetNumVerts();
2133 
2134  if(m_domainRange->m_doXrange)
2135  {
2136  int ncnt_low = 0;
2137  int ncnt_up = 0;
2138 
2139  for(int i = 0; i < nverts; ++i)
2140  {
2141  NekDouble xval = (*geom.GetVertex(i))[0];
2142  if(xval < m_domainRange->m_xmin)
2143  {
2144  ncnt_low++;
2145  }
2146 
2147  if(xval > m_domainRange->m_xmax)
2148  {
2149  ncnt_up++;
2150  }
2151  }
2152 
2153  // check for all verts to be less or greater than
2154  // range so that if element spans thin range then
2155  // it is still included
2156  if((ncnt_up == nverts)||(ncnt_low == nverts))
2157  {
2158  returnval = false;
2159  }
2160  }
2161 
2162  if(m_domainRange->m_doYrange)
2163  {
2164  int ncnt_low = 0;
2165  int ncnt_up = 0;
2166  for(int i = 0; i < nverts; ++i)
2167  {
2168  NekDouble yval = (*geom.GetVertex(i))[1];
2169  if(yval < m_domainRange->m_ymin)
2170  {
2171  ncnt_low++;
2172  }
2173 
2174  if(yval > m_domainRange->m_ymax)
2175  {
2176  ncnt_up++;
2177  }
2178  }
2179 
2180  // check for all verts to be less or greater than
2181  // range so that if element spans thin range then
2182  // it is still included
2183  if((ncnt_up == nverts)||(ncnt_low == nverts))
2184  {
2185  returnval = false;
2186  }
2187  }
2188 
2189  if(m_domainRange->m_doZrange)
2190  {
2191  int ncnt_low = 0;
2192  int ncnt_up = 0;
2193  for(int i = 0; i < nverts; ++i)
2194  {
2195  NekDouble zval = (*geom.GetVertex(i))[2];
2196 
2197  if(zval < m_domainRange->m_zmin)
2198  {
2199  ncnt_low++;
2200  }
2201 
2202  if(zval > m_domainRange->m_zmax)
2203  {
2204  ncnt_up++;
2205  }
2206  }
2207 
2208  // check for all verts to be less or greater than
2209  // range so that if element spans thin range then
2210  // it is still included
2211  if((ncnt_up == nverts)||(ncnt_low == nverts))
2212  {
2213  returnval = false;
2214  }
2215  }
2216 
2217  if(m_domainRange->m_checkShape)
2218  {
2219  if(geom.GetShapeType() != m_domainRange->m_shapeType)
2220  {
2221  returnval = false;
2222  }
2223  }
2224 
2225  }
2226 
2227  return returnval;
2228  }
2229 
2230  /**
2231  *
2232  */
2233  GeometrySharedPtr MeshGraph::GetCompositeItem(int whichComposite, int whichItem)
2234  {
2235  GeometrySharedPtr returnval;
2236  bool error = false;
2237 
2238  if (whichComposite >= 0 && whichComposite < int(m_meshComposites.size()))
2239  {
2240  if (whichItem >= 0 && whichItem < int(m_meshComposites[whichComposite]->size()))
2241  {
2242  returnval = m_meshComposites[whichComposite]->at(whichItem);
2243  }
2244  else
2245  {
2246  error = true;
2247  }
2248  }
2249  else
2250  {
2251  error = true;
2252  }
2253 
2254  if (error)
2255  {
2256  std::ostringstream errStream;
2257  errStream << "Unable to access composite item [" << whichComposite << "][" << whichItem << "].";
2258 
2259  std::string testStr = errStream.str();
2260 
2261  NEKERROR(ErrorUtil::efatal, testStr.c_str());
2262  }
2263 
2264  return returnval;
2265  }
2266 
2267 
2268  /**
2269  *
2270  */
2271  void MeshGraph::GetCompositeList(const std::string &compositeStr, CompositeMap &compositeVector) const
2272  {
2273  // Parse the composites into a list.
2274  typedef vector<unsigned int> SeqVector;
2275  SeqVector seqVector;
2276  bool parseGood = ParseUtils::GenerateSeqVector(compositeStr.c_str(), seqVector);
2277 
2278  ASSERTL0(parseGood && !seqVector.empty(), (std::string("Unable to read composite index range: ") + compositeStr).c_str());
2279 
2280  SeqVector addedVector; // Vector of those composites already added to compositeVector;
2281  for (SeqVector::iterator iter = seqVector.begin(); iter != seqVector.end(); ++iter)
2282  {
2283  // Only add a new one if it does not already exist in vector.
2284  // Can't go back and delete with a vector, so prevent it from
2285  // being added in the first place.
2286  if (std::find(addedVector.begin(), addedVector.end(), *iter) == addedVector.end())
2287  {
2288 
2289  // If the composite listed is not found and we are working
2290  // on a partitioned mesh, silently ignore it.
2291  if (m_meshComposites.find(*iter) == m_meshComposites.end()
2292  && m_meshPartitioned)
2293  {
2294  continue;
2295  }
2296 
2297  addedVector.push_back(*iter);
2298  Composite composite = GetComposite(*iter);
2299  CompositeMap::iterator compIter;
2300  if (composite)
2301  {
2302  compositeVector[*iter] = composite;
2303  }
2304  else
2305  {
2306  char str[64];
2307  ::sprintf(str, "%d", *iter);
2308  NEKERROR(ErrorUtil::ewarning, (std::string("Undefined composite: ") + str).c_str());
2309 
2310  }
2311  }
2312  }
2313  }
2314 
2315 
2316  /**
2317  *
2318  */
2319  const ExpansionMap &MeshGraph::GetExpansions(const std::string variable)
2320  {
2321  ExpansionMapShPtr returnval;
2322 
2323  if(m_expansionMapShPtrMap.count(variable))
2324  {
2325  returnval = m_expansionMapShPtrMap.find(variable)->second;
2326  }
2327  else
2328  {
2329  if(m_expansionMapShPtrMap.count("DefaultVar") == 0)
2330  {
2331  NEKERROR(ErrorUtil::efatal, (std::string("Unable to find expansion vector definition for field: ")+variable).c_str());
2332  }
2333  returnval = m_expansionMapShPtrMap.find("DefaultVar")->second;
2334  m_expansionMapShPtrMap[variable] = returnval;
2335 
2336  NEKERROR(ErrorUtil::ewarning, (std::string("Using Default variable expansion definition for field: ")+variable).c_str());
2337  }
2338 
2339  return *returnval;
2340  }
2341 
2342 
2343  /**
2344  *
2345  */
2346  ExpansionShPtr MeshGraph::GetExpansion(GeometrySharedPtr geom, const std::string variable)
2347  {
2348  ExpansionMapIter iter;
2349  ExpansionMapShPtr expansionMap = m_expansionMapShPtrMap.find(variable)->second;
2350 
2351  iter = expansionMap->find(geom->GetGlobalID());
2352  ASSERTL1(iter != expansionMap->end(),
2353  "Could not find expansion " +
2354  boost::lexical_cast<string>(geom->GetGlobalID()) +
2355  " in expansion for variable " + variable);
2356  return iter->second;
2357  }
2358 
2359 
2360  /**
2361  *
2362  */
2363  void MeshGraph::SetExpansions(
2364  std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
2365  {
2366  int i, j, k, cnt, id;
2367  GeometrySharedPtr geom;
2368 
2369  ExpansionMapShPtr expansionMap;
2370 
2371  // Loop over fields and determine unique fields string and
2372  // declare whole expansion list
2373  for(i = 0; i < fielddef.size(); ++i)
2374  {
2375  for(j = 0; j < fielddef[i]->m_fields.size(); ++j)
2376  {
2377  std::string field = fielddef[i]->m_fields[j];
2378  if(m_expansionMapShPtrMap.count(field) == 0)
2379  {
2380  expansionMap = SetUpExpansionMap();
2381  m_expansionMapShPtrMap[field] = expansionMap;
2382 
2383  // check to see if DefaultVar also not set and
2384  // if so assign it to this expansion
2385  if(m_expansionMapShPtrMap.count("DefaultVar") == 0)
2386  {
2387  m_expansionMapShPtrMap["DefaultVar"] = expansionMap;
2388  }
2389  }
2390  }
2391  }
2392 
2393  // loop over all elements find the geometry shared ptr and
2394  // set up basiskey vector
2395  for(i = 0; i < fielddef.size(); ++i)
2396  {
2397  cnt = 0;
2398  std::vector<std::string> fields = fielddef[i]->m_fields;
2399  std::vector<unsigned int> nmodes = fielddef[i]->m_numModes;
2400  std::vector<LibUtilities::BasisType> basis = fielddef[i]->m_basis;
2401  bool pointDef = fielddef[i]->m_pointsDef;
2402  bool numPointDef = fielddef[i]->m_numPointsDef;
2403 
2404  // Check points and numpoints
2405  std::vector<unsigned int> npoints = fielddef[i]->m_numPoints;
2406  std::vector<LibUtilities::PointsType> points = fielddef[i]->m_points;
2407 
2408  bool UniOrder = fielddef[i]->m_uniOrder;
2409 
2410  for (j = 0; j < fielddef[i]->m_elementIDs.size(); ++j)
2411  {
2412 
2413  LibUtilities::BasisKeyVector bkeyvec;
2414  id = fielddef[i]->m_elementIDs[j];
2415 
2416  switch (fielddef[i]->m_shapeType)
2417  {
2418  case LibUtilities::eSegment:
2419  {
2420  if(m_segGeoms.count(fielddef[i]->m_elementIDs[j]) == 0)
2421  {
2422  // skip element likely from parallel read
2423  if(!UniOrder)
2424  {
2425  cnt++;
2426  cnt += fielddef[i]->m_numHomogeneousDir;
2427  }
2428  continue;
2429  }
2430  geom = m_segGeoms[fielddef[i]->m_elementIDs[j]];
2431 
2432  LibUtilities::PointsKey pkey(nmodes[cnt]+1, LibUtilities::eGaussLobattoLegendre);
2433 
2434  if(numPointDef&&pointDef)
2435  {
2436  const LibUtilities::PointsKey pkey1(npoints[cnt], points[0]);
2437  pkey = pkey1;
2438  }
2439  else if(!numPointDef&&pointDef)
2440  {
2441  const LibUtilities::PointsKey pkey1(nmodes[cnt]+1, points[0]);
2442  pkey = pkey1;
2443  }
2444  else if(numPointDef&&!pointDef)
2445  {
2446  const LibUtilities::PointsKey pkey1(npoints[cnt], LibUtilities::eGaussLobattoLegendre);
2447  pkey = pkey1;
2448  }
2449 
2450  LibUtilities::BasisKey bkey(basis[0], nmodes[cnt], pkey);
2451 
2452  if(!UniOrder)
2453  {
2454  cnt++;
2455  cnt += fielddef[i]->m_numHomogeneousDir;
2456  }
2457  bkeyvec.push_back(bkey);
2458  }
2459  break;
2460  case LibUtilities::eTriangle:
2461  {
2462  if(m_triGeoms.count(fielddef[i]->m_elementIDs[j]) == 0)
2463  {
2464  // skip element likely from parallel read
2465  if(!UniOrder)
2466  {
2467  cnt += 2;
2468  cnt += fielddef[i]->m_numHomogeneousDir;
2469  }
2470  continue;
2471  }
2472  geom = m_triGeoms[fielddef[i]->m_elementIDs[j]];
2473 
2474  LibUtilities::PointsKey pkey(nmodes[cnt]+1, LibUtilities::eGaussLobattoLegendre);
2475  if(numPointDef&&pointDef)
2476  {
2477  const LibUtilities::PointsKey pkey2(npoints[cnt], points[0]);
2478  pkey = pkey2;
2479  }
2480  else if(!numPointDef&&pointDef)
2481  {
2482  const LibUtilities::PointsKey pkey2(nmodes[cnt]+1, points[0]);
2483  pkey = pkey2;
2484  }
2485  else if(numPointDef&&!pointDef)
2486  {
2487  const LibUtilities::PointsKey pkey2(npoints[cnt], LibUtilities::eGaussLobattoLegendre);
2488  pkey = pkey2;
2489  }
2490  LibUtilities::BasisKey bkey(basis[0], nmodes[cnt], pkey);
2491 
2492  bkeyvec.push_back(bkey);
2493 
2494  LibUtilities::PointsKey pkey1(nmodes[cnt+1], LibUtilities::eGaussRadauMAlpha1Beta0);
2495  if(numPointDef&&pointDef)
2496  {
2497  const LibUtilities::PointsKey pkey2(npoints[cnt+1], points[1]);
2498  pkey1 = pkey2;
2499  }
2500  else if(!numPointDef&&pointDef)
2501  {
2502  const LibUtilities::PointsKey pkey2(nmodes[cnt+1], points[1]);
2503  pkey1 = pkey2;
2504  }
2505  else if(numPointDef&&!pointDef)
2506  {
2507  const LibUtilities::PointsKey pkey2(npoints[cnt+1], LibUtilities::eGaussRadauMAlpha1Beta0);
2508  pkey1 = pkey2;
2509  }
2510  LibUtilities::BasisKey bkey1(basis[1], nmodes[cnt+1], pkey1);
2511  bkeyvec.push_back(bkey1);
2512 
2513  if(!UniOrder)
2514  {
2515  cnt += 2;
2516  cnt += fielddef[i]->m_numHomogeneousDir;
2517  }
2518  }
2519  break;
2520  case LibUtilities::eQuadrilateral:
2521  {
2522  if(m_quadGeoms.count(fielddef[i]->m_elementIDs[j]) == 0)
2523  {
2524  // skip element likely from parallel read
2525  if(!UniOrder)
2526  {
2527  cnt += 2;
2528  cnt += fielddef[i]->m_numHomogeneousDir;
2529  }
2530  continue;
2531  }
2532 
2533  geom = m_quadGeoms[fielddef[i]->m_elementIDs[j]];
2534 
2535  for(int b = 0; b < 2; ++b)
2536  {
2537  LibUtilities::PointsKey pkey(nmodes[cnt+b]+1, LibUtilities::eGaussLobattoLegendre);
2538 
2539  if(numPointDef&&pointDef)
2540  {
2541  const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
2542  pkey = pkey2;
2543  }
2544  else if(!numPointDef&&pointDef)
2545  {
2546  const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
2547  pkey = pkey2;
2548  }
2549  else if(numPointDef&&!pointDef)
2550  {
2551  const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
2552  pkey = pkey2;
2553  }
2554  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2555  bkeyvec.push_back(bkey);
2556  }
2557 
2558  if(!UniOrder)
2559  {
2560  cnt += 2;
2561  cnt += fielddef[i]->m_numHomogeneousDir;
2562  }
2563  }
2564  break;
2565 
2566  case LibUtilities::eTetrahedron:
2567  {
2568  k = fielddef[i]->m_elementIDs[j];
2569 
2570  // allow for possibility that fielddef is
2571  // larger than m_graph which can happen in
2572  // parallel runs
2573  if(m_tetGeoms.count(k) == 0)
2574  {
2575  if(!UniOrder)
2576  {
2577  cnt += 3;
2578  }
2579  continue;
2580  }
2581  geom = m_tetGeoms[k];
2582 
2583  {
2584  LibUtilities::PointsKey pkey(nmodes[cnt]+1, LibUtilities::eGaussLobattoLegendre);
2585 
2586  if(numPointDef&&pointDef)
2587  {
2588  const LibUtilities::PointsKey pkey2(npoints[cnt],points[0]);
2589  pkey = pkey2;
2590  }
2591  else if(!numPointDef&&pointDef)
2592  {
2593  const LibUtilities::PointsKey pkey2(nmodes[cnt]+1,points[0]);
2594  pkey = pkey2;
2595  }
2596  else if(numPointDef&&!pointDef)
2597  {
2598  const LibUtilities::PointsKey pkey2(npoints[cnt],LibUtilities::eGaussLobattoLegendre);
2599  pkey = pkey2;
2600  }
2601 
2602  LibUtilities::BasisKey bkey(basis[0],nmodes[cnt],pkey);
2603 
2604  bkeyvec.push_back(bkey);
2605  }
2606  {
2607  LibUtilities::PointsKey pkey(nmodes[cnt+1], LibUtilities::eGaussRadauMAlpha1Beta0);
2608 
2609  if(numPointDef&&pointDef)
2610  {
2611  const LibUtilities::PointsKey pkey2(npoints[cnt+1],points[1]);
2612  pkey = pkey2;
2613  }
2614  else if(!numPointDef&&pointDef)
2615  {
2616  const LibUtilities::PointsKey pkey2(nmodes[cnt+1]+1,points[1]);
2617  pkey = pkey2;
2618  }
2619  else if(numPointDef&&!pointDef)
2620  {
2621  const LibUtilities::PointsKey pkey2(npoints[cnt+1],LibUtilities::eGaussRadauMAlpha1Beta0);
2622  pkey = pkey2;
2623  }
2624 
2625  LibUtilities::BasisKey bkey(basis[1],nmodes[cnt+1],pkey);
2626 
2627  bkeyvec.push_back(bkey);
2628  }
2629 
2630  {
2631  LibUtilities::PointsKey pkey(nmodes[cnt+2], LibUtilities::eGaussRadauMAlpha2Beta0);
2632 
2633  if(numPointDef&&pointDef)
2634  {
2635  const LibUtilities::PointsKey pkey2(npoints[cnt+2],points[2]);
2636  pkey = pkey2;
2637  }
2638  else if(!numPointDef&&pointDef)
2639  {
2640  const LibUtilities::PointsKey pkey2(nmodes[cnt+2]+1,points[2]);
2641  pkey = pkey2;
2642  }
2643  else if(numPointDef&&!pointDef)
2644  {
2645  const LibUtilities::PointsKey pkey2(npoints[cnt+2],LibUtilities::eGaussRadauMAlpha1Beta0);
2646  pkey = pkey2;
2647  }
2648 
2649  LibUtilities::BasisKey bkey(basis[2],nmodes[cnt+2],pkey);
2650 
2651  bkeyvec.push_back(bkey);
2652  }
2653 
2654  if(!UniOrder)
2655  {
2656  cnt += 3;
2657  }
2658  }
2659  break;
2660  case LibUtilities::ePrism:
2661  {
2662  k = fielddef[i]->m_elementIDs[j];
2663  if(m_prismGeoms.count(k) == 0)
2664  {
2665  if(!UniOrder)
2666  {
2667  cnt += 3;
2668  }
2669  continue;
2670  }
2671  geom = m_prismGeoms[k];
2672 
2673  for(int b = 0; b < 2; ++b)
2674  {
2675  LibUtilities::PointsKey pkey(nmodes[cnt+b]+1,LibUtilities::eGaussLobattoLegendre);
2676 
2677  if(numPointDef&&pointDef)
2678  {
2679  const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
2680  pkey = pkey2;
2681  }
2682  else if(!numPointDef&&pointDef)
2683  {
2684  const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
2685  pkey = pkey2;
2686  }
2687  else if(numPointDef&&!pointDef)
2688  {
2689  const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
2690  pkey = pkey2;
2691  }
2692 
2693  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2694  bkeyvec.push_back(bkey);
2695  }
2696 
2697  {
2698  LibUtilities::PointsKey pkey(nmodes[cnt+2],LibUtilities::eGaussRadauMAlpha1Beta0);
2699 
2700  if(numPointDef&&pointDef)
2701  {
2702  const LibUtilities::PointsKey pkey2(npoints[cnt+2],points[2]);
2703  pkey = pkey2;
2704  }
2705  else if(!numPointDef&&pointDef)
2706  {
2707  const LibUtilities::PointsKey pkey2(nmodes[cnt+2]+1,points[2]);
2708  pkey = pkey2;
2709  }
2710  else if(numPointDef&&!pointDef)
2711  {
2712  const LibUtilities::PointsKey pkey2(npoints[cnt+2],LibUtilities::eGaussLobattoLegendre);
2713  pkey = pkey2;
2714  }
2715 
2716  LibUtilities::BasisKey bkey(basis[2],nmodes[cnt+2],pkey);
2717  bkeyvec.push_back(bkey);
2718  }
2719 
2720  if(!UniOrder)
2721  {
2722  cnt += 3;
2723  }
2724  }
2725  break;
2726  case LibUtilities::ePyramid:
2727  {
2728  k = fielddef[i]->m_elementIDs[j];
2729  ASSERTL0(m_pyrGeoms.find(k) != m_pyrGeoms.end(),
2730  "Failed to find geometry with same global id");
2731  geom = m_pyrGeoms[k];
2732 
2733 
2734  for(int b = 0; b < 2; ++b)
2735  {
2736  LibUtilities::PointsKey pkey(nmodes[cnt+b]+1,LibUtilities::eGaussLobattoLegendre);
2737 
2738  if(numPointDef&&pointDef)
2739  {
2740  const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
2741  pkey = pkey2;
2742  }
2743  else if(!numPointDef&&pointDef)
2744  {
2745  const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
2746  pkey = pkey2;
2747  }
2748  else if(numPointDef&&!pointDef)
2749  {
2750  const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
2751  pkey = pkey2;
2752  }
2753 
2754  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2755  bkeyvec.push_back(bkey);
2756  }
2757 
2758  {
2759  LibUtilities::PointsKey pkey(nmodes[cnt+2],LibUtilities::eGaussRadauMAlpha2Beta0);
2760 
2761  if(numPointDef&&pointDef)
2762  {
2763  const LibUtilities::PointsKey pkey2(npoints[cnt+2],points[2]);
2764  pkey = pkey2;
2765  }
2766  else if(!numPointDef&&pointDef)
2767  {
2768  const LibUtilities::PointsKey pkey2(nmodes[cnt+2]+1,points[2]);
2769  pkey = pkey2;
2770  }
2771  else if(numPointDef&&!pointDef)
2772  {
2773  const LibUtilities::PointsKey pkey2(npoints[cnt+2],LibUtilities::eGaussLobattoLegendre);
2774  pkey = pkey2;
2775  }
2776 
2777  LibUtilities::BasisKey bkey(basis[2],nmodes[cnt+2],pkey);
2778  bkeyvec.push_back(bkey);
2779  }
2780 
2781  if(!UniOrder)
2782  {
2783  cnt += 3;
2784  }
2785  }
2786  break;
2787  case LibUtilities::eHexahedron:
2788  {
2789  k = fielddef[i]->m_elementIDs[j];
2790  if(m_hexGeoms.count(k) == 0)
2791  {
2792  if(!UniOrder)
2793  {
2794  cnt += 3;
2795  }
2796  continue;
2797  }
2798 
2799  geom = m_hexGeoms[k];
2800 
2801  for(int b = 0; b < 3; ++b)
2802  {
2803  LibUtilities::PointsKey pkey(nmodes[cnt+b],LibUtilities::eGaussLobattoLegendre);
2804 
2805  if(numPointDef&&pointDef)
2806  {
2807  const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
2808  pkey = pkey2;
2809  }
2810  else if(!numPointDef&&pointDef)
2811  {
2812  const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
2813  pkey = pkey2;
2814  }
2815  else if(numPointDef&&!pointDef)
2816  {
2817  const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
2818  pkey = pkey2;
2819  }
2820 
2821  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2822  bkeyvec.push_back(bkey);
2823  }
2824 
2825  if(!UniOrder)
2826  {
2827  cnt += 3;
2828  }
2829  }
2830  break;
2831  default:
2832  ASSERTL0(false,"Need to set up for pyramid and prism 3D Expansions");
2833  break;
2834  }
2835 
2836  for(k = 0; k < fields.size(); ++k)
2837  {
2838  expansionMap = m_expansionMapShPtrMap.find(fields[k])->second;
2839  if((*expansionMap).find(id) != (*expansionMap).end())
2840  {
2841  (*expansionMap)[id]->m_geomShPtr = geom;
2842  (*expansionMap)[id]->m_basisKeyVector = bkeyvec;
2843  }
2844  }
2845  }
2846  }
2847  }
2848 
2849 
2850  /**
2851  *
2852  */
2853  void MeshGraph::SetExpansions(
2854  std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef,
2855  std::vector< std::vector<LibUtilities::PointsType> > &pointstype)
2856  {
2857  int i,j,k,cnt,id;
2858  GeometrySharedPtr geom;
2859 
2860  ExpansionMapShPtr expansionMap;
2861 
2862  // Loop over fields and determine unique fields string and
2863  // declare whole expansion list
2864  for(i = 0; i < fielddef.size(); ++i)
2865  {
2866  for(j = 0; j < fielddef[i]->m_fields.size(); ++j)
2867  {
2868  std::string field = fielddef[i]->m_fields[j];
2869  if(m_expansionMapShPtrMap.count(field) == 0)
2870  {
2871  expansionMap = SetUpExpansionMap();
2872  m_expansionMapShPtrMap[field] = expansionMap;
2873 
2874  // check to see if DefaultVar also not set and
2875  // if so assign it to this expansion
2876  if(m_expansionMapShPtrMap.count("DefaultVar") == 0)
2877  {
2878  m_expansionMapShPtrMap["DefaultVar"] = expansionMap;
2879  }
2880  }
2881  }
2882  }
2883 
2884 
2885  // loop over all elements find the geometry shared ptr and
2886  // set up basiskey vector
2887  for(i = 0; i < fielddef.size(); ++i)
2888  {
2889  cnt = 0;
2890  std::vector<std::string> fields = fielddef[i]->m_fields;
2891  std::vector<unsigned int> nmodes = fielddef[i]->m_numModes;
2892  std::vector<LibUtilities::BasisType> basis = fielddef[i]->m_basis;
2893  bool UniOrder = fielddef[i]->m_uniOrder;
2894 
2895  for(j = 0; j < fielddef[i]->m_elementIDs.size(); ++j)
2896  {
2897  LibUtilities::BasisKeyVector bkeyvec;
2898  id = fielddef[i]->m_elementIDs[j];
2899 
2900  switch(fielddef[i]->m_shapeType)
2901  {
2902  case LibUtilities::eSegment:
2903  {
2904  k = fielddef[i]->m_elementIDs[j];
2905  ASSERTL0(m_segGeoms.find(k) != m_segGeoms.end(),
2906  "Failed to find geometry with same global id.");
2907  geom = m_segGeoms[k];
2908 
2909  const LibUtilities::PointsKey pkey(nmodes[cnt], pointstype[i][0]);
2910  LibUtilities::BasisKey bkey(basis[0], nmodes[cnt], pkey);
2911  if(!UniOrder)
2912  {
2913  cnt++;
2914  cnt += fielddef[i]->m_numHomogeneousDir;
2915  }
2916  bkeyvec.push_back(bkey);
2917  }
2918  break;
2919  case LibUtilities::eTriangle:
2920  {
2921  k = fielddef[i]->m_elementIDs[j];
2922  ASSERTL0(m_triGeoms.find(k) != m_triGeoms.end(),
2923  "Failed to find geometry with same global id.");
2924  geom = m_triGeoms[k];
2925  for(int b = 0; b < 2; ++b)
2926  {
2927  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
2928  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2929  bkeyvec.push_back(bkey);
2930  }
2931 
2932  if(!UniOrder)
2933  {
2934  cnt += 2;
2935  cnt += fielddef[i]->m_numHomogeneousDir;
2936  }
2937  }
2938  break;
2939  case LibUtilities::eQuadrilateral:
2940  {
2941  k = fielddef[i]->m_elementIDs[j];
2942  ASSERTL0(m_quadGeoms.find(k) != m_quadGeoms.end(),
2943  "Failed to find geometry with same global id");
2944  geom = m_quadGeoms[k];
2945 
2946  for(int b = 0; b < 2; ++b)
2947  {
2948  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
2949  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2950  bkeyvec.push_back(bkey);
2951  }
2952 
2953  if(!UniOrder)
2954  {
2955  cnt += 2;
2956  cnt += fielddef[i]->m_numHomogeneousDir;
2957  }
2958  }
2959  break;
2960  case LibUtilities::eTetrahedron:
2961  {
2962  k = fielddef[i]->m_elementIDs[j];
2963  ASSERTL0(m_tetGeoms.find(k) != m_tetGeoms.end(),
2964  "Failed to find geometry with same global id");
2965  geom = m_tetGeoms[k];
2966 
2967  for(int b = 0; b < 3; ++b)
2968  {
2969  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
2970  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2971  bkeyvec.push_back(bkey);
2972  }
2973 
2974  if(!UniOrder)
2975  {
2976  cnt += 3;
2977  }
2978  }
2979  break;
2980  case LibUtilities::ePyramid:
2981  {
2982  k = fielddef[i]->m_elementIDs[j];
2983  ASSERTL0(m_pyrGeoms.find(k) != m_pyrGeoms.end(),
2984  "Failed to find geometry with same global id");
2985  geom = m_pyrGeoms[k];
2986 
2987  for(int b = 0; b < 3; ++b)
2988  {
2989  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
2990  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2991  bkeyvec.push_back(bkey);
2992  }
2993 
2994  if(!UniOrder)
2995  {
2996  cnt += 3;
2997  }
2998  }
2999  break;
3000  case LibUtilities::ePrism:
3001  {
3002  k = fielddef[i]->m_elementIDs[j];
3003  ASSERTL0(m_prismGeoms.find(k) != m_prismGeoms.end(),
3004  "Failed to find geometry with same global id");
3005  geom = m_prismGeoms[k];
3006 
3007  for(int b = 0; b < 3; ++b)
3008  {
3009  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
3010  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
3011  bkeyvec.push_back(bkey);
3012  }
3013 
3014  if(!UniOrder)
3015  {
3016  cnt += 3;
3017  }
3018  }
3019  break;
3020  case LibUtilities::eHexahedron:
3021  {
3022  k = fielddef[i]->m_elementIDs[j];
3023  ASSERTL0(m_hexGeoms.find(k) != m_hexGeoms.end(),
3024  "Failed to find geometry with same global id");
3025  geom = m_hexGeoms[k];
3026 
3027  for(int b = 0; b < 3; ++b)
3028  {
3029  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
3030  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
3031  bkeyvec.push_back(bkey);
3032  }
3033 
3034  if(!UniOrder)
3035  {
3036  cnt += 3;
3037  }
3038  }
3039  break;
3040  default:
3041  ASSERTL0(false,"Need to set up for pyramid and prism 3D Expansions");
3042  break;
3043  }
3044 
3045  for(k = 0; k < fields.size(); ++k)
3046  {
3047  expansionMap = m_expansionMapShPtrMap.find(fields[k])->second;
3048  if((*expansionMap).find(id) != (*expansionMap).end())
3049  {
3050  (*expansionMap)[id]->m_geomShPtr = geom;
3051  (*expansionMap)[id]->m_basisKeyVector = bkeyvec;
3052  }
3053  }
3054  }
3055  }
3056  }
3057 
3058  /**
3059  * \brief Reset all points keys to have equispaced points with
3060  * optional arguemt of \a npoints which redefines how many
3061  * points are to be used.
3062  */
3063  void MeshGraph::SetExpansionsToEvenlySpacedPoints(int npoints)
3064  {
3065  ExpansionMapShPtrMapIter it;
3066 
3067  // iterate over all defined expansions
3068  for(it = m_expansionMapShPtrMap.begin(); it != m_expansionMapShPtrMap.end(); ++it)
3069  {
3070  ExpansionMapIter expIt;
3071 
3072  for(expIt = it->second->begin(); expIt != it->second->end(); ++expIt)
3073  {
3074  for(int i = 0; i < expIt->second->m_basisKeyVector.size(); ++i)
3075  {
3076  LibUtilities::BasisKey bkeyold = expIt->second->m_basisKeyVector[i];
3077 
3078  int npts;
3079 
3080  if(npoints) // use input
3081  {
3082  npts = npoints;
3083  }
3084  else
3085  {
3086  npts = bkeyold.GetNumModes();
3087  }
3088 
3089 
3090  const LibUtilities::PointsKey pkey(npts,LibUtilities::ePolyEvenlySpaced);
3091  LibUtilities::BasisKey bkeynew(bkeyold.GetBasisType(),bkeyold.GetNumModes(), pkey);
3092  expIt->second->m_basisKeyVector[i] = bkeynew;
3093 
3094  }
3095  }
3096  }
3097  }
3098 
3099  /**
3100  * \brief Reset all points keys to have expansion order of \a
3101  * nmodes. we keep the point distribution the same and make
3102  * the number of points the same difference from the number
3103  * of modes as the original expansion definition
3104  */
3105  void MeshGraph::SetExpansionsToPolyOrder(int nmodes)
3106  {
3107  ExpansionMapShPtrMapIter it;
3108 
3109  // iterate over all defined expansions
3110  for(it = m_expansionMapShPtrMap.begin(); it != m_expansionMapShPtrMap.end(); ++it)
3111  {
3112  ExpansionMapIter expIt;
3113 
3114  for(expIt = it->second->begin(); expIt != it->second->end(); ++expIt)
3115  {
3116  for(int i = 0; i < expIt->second->m_basisKeyVector.size(); ++i)
3117  {
3118  LibUtilities::BasisKey bkeyold = expIt->second->m_basisKeyVector[i];
3119 
3120  int npts = nmodes + (bkeyold.GetNumPoints() - bkeyold.GetNumModes());
3121 
3122  const LibUtilities::PointsKey pkey(npts,bkeyold.GetPointsType());
3123  LibUtilities::BasisKey bkeynew(bkeyold.GetBasisType(),nmodes, pkey);
3124  expIt->second->m_basisKeyVector[i] = bkeynew;
3125 
3126  }
3127  }
3128  }
3129  }
3130 
3131 
3132  /**
3133  * \brief Reset all points keys to have expansion order of \a
3134  * nmodes. we keep the point distribution the same and make
3135  * the number of points the same difference from the number
3136  * of modes as the original expansion definition
3137  */
3138  void MeshGraph::SetExpansionsToPointOrder(int npts)
3139  {
3140  ExpansionMapShPtrMapIter it;
3141 
3142  // iterate over all defined expansions
3143  for (it = m_expansionMapShPtrMap.begin();
3144  it != m_expansionMapShPtrMap.end();
3145  ++it)
3146  {
3147  ExpansionMapIter expIt;
3148 
3149  for (expIt = it->second->begin();
3150  expIt != it->second->end();
3151  ++expIt)
3152  {
3153  for(int i = 0;
3154  i < expIt->second->m_basisKeyVector.size();
3155  ++i)
3156  {
3157  LibUtilities::BasisKey bkeyold =
3158  expIt->second->m_basisKeyVector[i];
3159 
3160  const LibUtilities::PointsKey pkey(
3161  npts, bkeyold.GetPointsType());
3162 
3163  LibUtilities::BasisKey bkeynew(bkeyold.GetBasisType(),
3164  bkeyold.GetNumModes(),
3165  pkey);
3166  expIt->second->m_basisKeyVector[i] = bkeynew;
3167  }
3168  }
3169  }
3170  }
3171 
3172 
3173  /**
3174  * For each element of shape given by \a shape in field \a
3175  * var, replace the current BasisKeyVector describing the
3176  * expansion in each dimension, with the one provided by \a
3177  * keys.
3178  *
3179  * @TODO: Allow selection of elements through a CompositeVector,
3180  * as well as by type.
3181  *
3182  * @param shape The shape of elements to be changed.
3183  * @param keys The new basis vector to apply to those elements.
3184  */
3185  void MeshGraph::SetBasisKey(LibUtilities::ShapeType shape,
3186  LibUtilities::BasisKeyVector &keys,
3187  std::string var)
3188  {
3189  ExpansionMapIter elemIter;
3190 
3191  ExpansionMapShPtr expansionMap = m_expansionMapShPtrMap.find(var)->second;
3192 
3193  for (elemIter = expansionMap->begin(); elemIter != expansionMap->end(); ++elemIter)
3194  {
3195  if ((elemIter->second)->m_geomShPtr->GetShapeType() == shape)
3196  {
3197  (elemIter->second)->m_basisKeyVector = keys;
3198  }
3199  }
3200  }
3201 
3202 
3203  /**
3204  *
3205  */
3206  LibUtilities::BasisKeyVector MeshGraph::DefineBasisKeyFromExpansionType(
3207  GeometrySharedPtr in,
3208  ExpansionType type,
3209  const int nummodes)
3210  {
3211  LibUtilities::BasisKeyVector returnval;
3212 
3213  LibUtilities::ShapeType shape= in->GetShapeType();
3214 
3215  int quadoffset = 1;
3216  switch(type)
3217  {
3218  case eModified:
3219  case eModifiedGLLRadau10:
3220  quadoffset = 1;
3221  break;
3222  case eModifiedQuadPlus1:
3223  quadoffset = 2;
3224  break;
3225  case eModifiedQuadPlus2:
3226  quadoffset = 3;
3227  break;
3228  default:
3229  break;
3230  }
3231 
3232  switch(type)
3233  {
3234  case eModified:
3235  case eModifiedQuadPlus1:
3236  case eModifiedQuadPlus2:
3237  case eModifiedGLLRadau10:
3238  {
3239  switch (shape)
3240  {
3241  case LibUtilities::eSegment:
3242  {
3243  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3244  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3245  returnval.push_back(bkey);
3246  }
3247  break;
3248  case LibUtilities::eQuadrilateral:
3249  {
3250  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3251  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3252  returnval.push_back(bkey);
3253  returnval.push_back(bkey);
3254  }
3255  break;
3256  case LibUtilities::eHexahedron:
3257  {
3258  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3259  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3260  returnval.push_back(bkey);
3261  returnval.push_back(bkey);
3262  returnval.push_back(bkey);
3263  }
3264  break;
3265  case LibUtilities::eTriangle:
3266  {
3267  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3268  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3269  returnval.push_back(bkey);
3270 
3271  const LibUtilities::PointsKey pkey1(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha1Beta0);
3272  LibUtilities::BasisKey bkey1(LibUtilities::eModified_B, nummodes, pkey1);
3273 
3274  returnval.push_back(bkey1);
3275  }
3276  break;
3277  case LibUtilities::eTetrahedron:
3278  {
3279  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3280  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3281  returnval.push_back(bkey);
3282 
3283  const LibUtilities::PointsKey pkey1(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha1Beta0);
3284  LibUtilities::BasisKey bkey1(LibUtilities::eModified_B, nummodes, pkey1);
3285  returnval.push_back(bkey1);
3286 
3287  if(type == eModifiedGLLRadau10)
3288  {
3289  const LibUtilities::PointsKey pkey2(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha1Beta0);
3290  LibUtilities::BasisKey bkey2(LibUtilities::eModified_C, nummodes, pkey2);
3291  returnval.push_back(bkey2);
3292  }
3293  else
3294  {
3295  const LibUtilities::PointsKey pkey2(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha2Beta0);
3296  LibUtilities::BasisKey bkey2(LibUtilities::eModified_C, nummodes, pkey2);
3297  returnval.push_back(bkey2);
3298  }
3299  }
3300  break;
3301  case LibUtilities::ePyramid:
3302  {
3303  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3304  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3305  returnval.push_back(bkey);
3306  returnval.push_back(bkey);
3307 
3308  const LibUtilities::PointsKey pkey1(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha2Beta0);
3309  LibUtilities::BasisKey bkey1(LibUtilities::eModified_C, nummodes, pkey1);
3310  returnval.push_back(bkey1);
3311  }
3312  break;
3313  case LibUtilities::ePrism:
3314  {
3315  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3316  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3317  returnval.push_back(bkey);
3318  returnval.push_back(bkey);
3319 
3320  const LibUtilities::PointsKey pkey1(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha1Beta0);
3321  LibUtilities::BasisKey bkey1(LibUtilities::eModified_B, nummodes, pkey1);
3322  returnval.push_back(bkey1);
3323 
3324  }
3325  break;
3326  default:
3327  {
3328  ASSERTL0(false,"Expansion not defined in switch for this shape");
3329  }
3330  break;
3331  }
3332  }
3333  break;
3334 
3335  case eGLL_Lagrange:
3336  {
3337  switch(shape)
3338  {
3339  case LibUtilities::eSegment:
3340  {
3341  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3342  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3343  returnval.push_back(bkey);
3344  }
3345  break;
3346  case LibUtilities::eQuadrilateral:
3347  {
3348  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3349  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3350  returnval.push_back(bkey);
3351  returnval.push_back(bkey);
3352  }
3353  break;
3354  case LibUtilities::eTriangle: // define with corrects points key
3355  // and change to Ortho on construction
3356  {
3357  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3358  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3359  returnval.push_back(bkey);
3360 
3361  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussRadauMAlpha1Beta0);
3362  LibUtilities::BasisKey bkey1(LibUtilities::eOrtho_B, nummodes, pkey1);
3363  returnval.push_back(bkey1);
3364  }
3365  break;
3366  case LibUtilities::eHexahedron:
3367  {
3368  const LibUtilities::PointsKey pkey(nummodes+1,LibUtilities::eGaussLobattoLegendre);
3369  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3370 
3371  returnval.push_back(bkey);
3372  returnval.push_back(bkey);
3373  returnval.push_back(bkey);
3374  }
3375  break;
3376  default:
3377  {
3378  ASSERTL0(false, "Expansion not defined in switch for this shape");
3379  }
3380  break;
3381  }
3382  }
3383  break;
3384 
3385  case eGauss_Lagrange:
3386  {
3387  switch (shape)
3388  {
3389  case LibUtilities::eSegment:
3390  {
3391  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussLegendre);
3392  LibUtilities::BasisKey bkey(LibUtilities::eGauss_Lagrange, nummodes, pkey);
3393 
3394  returnval.push_back(bkey);
3395  }
3396  break;
3397  case LibUtilities::eQuadrilateral:
3398  {
3399  const LibUtilities::PointsKey pkey(nummodes,LibUtilities::eGaussGaussLegendre);
3400  LibUtilities::BasisKey bkey(LibUtilities::eGauss_Lagrange, nummodes, pkey);
3401 
3402  returnval.push_back(bkey);
3403  returnval.push_back(bkey);
3404  }
3405  break;
3406  case LibUtilities::eHexahedron:
3407  {
3408  const LibUtilities::PointsKey pkey(nummodes,LibUtilities::eGaussGaussLegendre);
3409  LibUtilities::BasisKey bkey(LibUtilities::eGauss_Lagrange, nummodes, pkey);
3410 
3411  returnval.push_back(bkey);
3412  returnval.push_back(bkey);
3413  returnval.push_back(bkey);
3414  }
3415  break;
3416  default:
3417  {
3418  ASSERTL0(false, "Expansion not defined in switch for this shape");
3419  }
3420  break;
3421  }
3422  }
3423  break;
3424 
3425  case eOrthogonal:
3426  {
3427  switch (shape)
3428  {
3429  case LibUtilities::eSegment:
3430  {
3431  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3432  LibUtilities::BasisKey bkey(LibUtilities::eOrtho_A, nummodes, pkey);
3433 
3434  returnval.push_back(bkey);
3435  }
3436  break;
3437  case LibUtilities::eTriangle:
3438  {
3439  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3440  LibUtilities::BasisKey bkey(LibUtilities::eOrtho_A, nummodes, pkey);
3441 
3442  returnval.push_back(bkey);
3443 
3444  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussRadauMAlpha1Beta0);
3445  LibUtilities::BasisKey bkey1(LibUtilities::eOrtho_B, nummodes, pkey1);
3446 
3447  returnval.push_back(bkey1);
3448  }
3449  break;
3450  case LibUtilities::eQuadrilateral:
3451  {
3452  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3453  LibUtilities::BasisKey bkey(LibUtilities::eOrtho_A, nummodes, pkey);
3454 
3455  returnval.push_back(bkey);
3456  returnval.push_back(bkey);
3457  }
3458  break;
3459  case LibUtilities::eTetrahedron:
3460  {
3461  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3462  LibUtilities::BasisKey bkey(LibUtilities::eOrtho_A, nummodes, pkey);
3463 
3464  returnval.push_back(bkey);
3465 
3466  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussRadauMAlpha1Beta0);
3467  LibUtilities::BasisKey bkey1(LibUtilities::eOrtho_B, nummodes, pkey1);
3468 
3469  returnval.push_back(bkey1);
3470 
3471  const LibUtilities::PointsKey pkey2(nummodes, LibUtilities::eGaussRadauMAlpha2Beta0);
3472  LibUtilities::BasisKey bkey2(LibUtilities::eOrtho_C, nummodes, pkey2);
3473  }
3474  break;
3475  default:
3476  {
3477  ASSERTL0(false,"Expansion not defined in switch for this shape");
3478  }
3479  break;
3480  }
3481  }
3482  break;
3483 
3484  case eGLL_Lagrange_SEM:
3485  {
3486  switch (shape)
3487  {
3488  case LibUtilities::eSegment:
3489  {
3490  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussLobattoLegendre);
3491  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3492 
3493  returnval.push_back(bkey);
3494  }
3495  break;
3496  case LibUtilities::eQuadrilateral:
3497  {
3498  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussLobattoLegendre);
3499  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3500 
3501  returnval.push_back(bkey);
3502  returnval.push_back(bkey);
3503  }
3504  break;
3505  case LibUtilities::eHexahedron:
3506  {
3507  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussLobattoLegendre);
3508  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3509 
3510  returnval.push_back(bkey);
3511  returnval.push_back(bkey);
3512  returnval.push_back(bkey);
3513  }
3514  break;
3515  default:
3516  {
3517  ASSERTL0(false,"Expansion not defined in switch for this shape");
3518  }
3519  break;
3520  }
3521  }
3522  break;
3523 
3524 
3525  case eFourier:
3526  {
3527  switch (shape)
3528  {
3529  case LibUtilities::eSegment:
3530  {
3531  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3532  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3533  returnval.push_back(bkey);
3534  }
3535  break;
3536  case LibUtilities::eQuadrilateral:
3537  {
3538  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3539  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3540  returnval.push_back(bkey);
3541  returnval.push_back(bkey);
3542  }
3543  break;
3544  case LibUtilities::eHexahedron:
3545  {
3546  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3547  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3548  returnval.push_back(bkey);
3549  returnval.push_back(bkey);
3550  returnval.push_back(bkey);
3551  }
3552  break;
3553  default:
3554  {
3555  ASSERTL0(false,"Expansion not defined in switch for this shape");
3556  }
3557  break;
3558  }
3559  }
3560  break;
3561 
3562 
3563  case eFourierSingleMode:
3564  {
3565  switch (shape)
3566  {
3567  case LibUtilities::eSegment:
3568  {
3569  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3570  LibUtilities::BasisKey bkey(LibUtilities::eFourierSingleMode, nummodes, pkey);
3571  returnval.push_back(bkey);
3572  }
3573  break;
3574  case LibUtilities::eQuadrilateral:
3575  {
3576  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3577  LibUtilities::BasisKey bkey(LibUtilities::eFourierSingleMode, nummodes, pkey);
3578  returnval.push_back(bkey);
3579  returnval.push_back(bkey);
3580  }
3581  break;
3582  case LibUtilities::eHexahedron:
3583  {
3584  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3585  LibUtilities::BasisKey bkey(LibUtilities::eFourierSingleMode, nummodes, pkey);
3586  returnval.push_back(bkey);
3587  returnval.push_back(bkey);
3588  returnval.push_back(bkey);
3589  }
3590  break;
3591  default:
3592  {
3593  ASSERTL0(false,"Expansion not defined in switch for this shape");
3594  }
3595  break;
3596  }
3597  }
3598  break;
3599 
3600  case eFourierHalfModeRe:
3601  {
3602  switch (shape)
3603  {
3604  case LibUtilities::eSegment:
3605  {
3606  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3607  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeRe, nummodes, pkey);
3608  returnval.push_back(bkey);
3609  }
3610  break;
3611  case LibUtilities::eQuadrilateral:
3612  {
3613  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3614  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeRe, nummodes, pkey);
3615  returnval.push_back(bkey);
3616  returnval.push_back(bkey);
3617  }
3618  break;
3619  case LibUtilities::eHexahedron:
3620  {
3621  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3622  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeRe, nummodes, pkey);
3623  returnval.push_back(bkey);
3624  returnval.push_back(bkey);
3625  returnval.push_back(bkey);
3626  }
3627  break;
3628  default:
3629  {
3630  ASSERTL0(false,"Expansion not defined in switch for this shape");
3631  }
3632  break;
3633  }
3634  }
3635  break;
3636 
3637  case eFourierHalfModeIm:
3638  {
3639  switch (shape)
3640  {
3641  case LibUtilities::eSegment:
3642  {
3643  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3644  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeIm, nummodes, pkey);
3645  returnval.push_back(bkey);
3646  }
3647  break;
3648  case LibUtilities::eQuadrilateral:
3649  {
3650  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3651  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeIm, nummodes, pkey);
3652  returnval.push_back(bkey);
3653  returnval.push_back(bkey);
3654  }
3655  break;
3656  case LibUtilities::eHexahedron:
3657  {
3658  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3659  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeIm, nummodes, pkey);
3660  returnval.push_back(bkey);
3661  returnval.push_back(bkey);
3662  returnval.push_back(bkey);
3663  }
3664  break;
3665  default:
3666  {
3667  ASSERTL0(false,"Expansion not defined in switch for this shape");
3668  }
3669  break;
3670  }
3671  }
3672  break;
3673 
3674  case eChebyshev:
3675  {
3676  switch (shape)
3677  {
3678  case LibUtilities::eSegment:
3679  {
3680  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussChebyshev);
3681  LibUtilities::BasisKey bkey(LibUtilities::eChebyshev, nummodes, pkey);
3682  returnval.push_back(bkey);
3683  }
3684  break;
3685  case LibUtilities::eQuadrilateral:
3686  {
3687  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussChebyshev);
3688  LibUtilities::BasisKey bkey(LibUtilities::eChebyshev, nummodes, pkey);
3689  returnval.push_back(bkey);
3690  returnval.push_back(bkey);
3691  }
3692  break;
3693  case LibUtilities::eHexahedron:
3694  {
3695  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussChebyshev);
3696  LibUtilities::BasisKey bkey(LibUtilities::eChebyshev, nummodes, pkey);
3697  returnval.push_back(bkey);
3698  returnval.push_back(bkey);
3699  returnval.push_back(bkey);
3700  }
3701  break;
3702  default:
3703  {
3704  ASSERTL0(false,"Expansion not defined in switch for this shape");
3705  }
3706  break;
3707  }
3708  }
3709  break;
3710 
3711  case eFourierChebyshev:
3712  {
3713  switch (shape)
3714  {
3715  case LibUtilities::eQuadrilateral:
3716  {
3717  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3718  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3719  returnval.push_back(bkey);
3720 
3721  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussGaussChebyshev);
3722  LibUtilities::BasisKey bkey1(LibUtilities::eChebyshev, nummodes, pkey1);
3723  returnval.push_back(bkey1);
3724  }
3725  break;
3726  default:
3727  {
3728  ASSERTL0(false,"Expansion not defined in switch for this shape");
3729  }
3730  break;
3731  }
3732  }
3733  break;
3734 
3735  case eChebyshevFourier:
3736  {
3737  switch (shape)
3738  {
3739  case LibUtilities::eQuadrilateral:
3740  {
3741  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussGaussChebyshev);
3742  LibUtilities::BasisKey bkey1(LibUtilities::eChebyshev, nummodes, pkey1);
3743  returnval.push_back(bkey1);
3744 
3745  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3746  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3747  returnval.push_back(bkey);
3748  }
3749  break;
3750  default:
3751  {
3752  ASSERTL0(false,"Expansion not defined in switch for this shape");
3753  }
3754  break;
3755  }
3756  }
3757  break;
3758 
3759  case eFourierModified:
3760  {
3761  switch (shape)
3762  {
3763  case LibUtilities::eQuadrilateral:
3764  {
3765  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3766  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3767  returnval.push_back(bkey);
3768 
3769  const LibUtilities::PointsKey pkey1(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3770  LibUtilities::BasisKey bkey1(LibUtilities::eModified_A, nummodes, pkey1);
3771  returnval.push_back(bkey1);
3772  }
3773  break;
3774  default:
3775  {
3776  ASSERTL0(false,"Expansion not defined in switch for this shape");
3777  }
3778  break;
3779  }
3780  }
3781  break;
3782 
3783  default:
3784  {
3785  ASSERTL0(false,"Expansion type not defined");
3786  }
3787  break;
3788 
3789  }
3790 
3791  return returnval;
3792  }
3793 
3794  /**
3795  *
3796  */
3797  LibUtilities::BasisKeyVector
3798  MeshGraph::DefineBasisKeyFromExpansionTypeHomo(
3799  GeometrySharedPtr in,
3800  ExpansionType type_x,
3801  ExpansionType type_y,
3802  ExpansionType type_z,
3803  const int nummodes_x,
3804  const int nummodes_y,
3805  const int nummodes_z)
3806  {
3807  LibUtilities::BasisKeyVector returnval;
3808 
3809  LibUtilities::ShapeType shape = in->GetShapeType();
3810 
3811  switch (shape)
3812  {
3813  case LibUtilities::eSegment:
3814  {
3815  ASSERTL0(false,"Homogeneous expansion not defined for this shape");
3816  }
3817  break;
3818 
3819  case LibUtilities::eQuadrilateral:
3820  {
3821  ASSERTL0(false,"Homogeneous expansion not defined for this shape");
3822  }
3823  break;
3824 
3825  case LibUtilities::eHexahedron:
3826  {
3827  switch(type_x)
3828  {
3829  case eFourier:
3830  {
3831  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eFourierEvenlySpaced);
3832  LibUtilities::BasisKey bkey1(LibUtilities::eFourier,nummodes_x,pkey1);
3833  returnval.push_back(bkey1);
3834  }
3835  break;
3836 
3837  case eFourierSingleMode:
3838  {
3839  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eFourierSingleModeSpaced);
3840  LibUtilities::BasisKey bkey1(LibUtilities::eFourierSingleMode,nummodes_x,pkey1);
3841  returnval.push_back(bkey1);
3842  }
3843  break;
3844 
3845  case eFourierHalfModeRe:
3846  {
3847  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eFourierSingleModeSpaced);
3848  LibUtilities::BasisKey bkey1(LibUtilities::eFourierHalfModeRe,nummodes_x,pkey1);
3849  returnval.push_back(bkey1);
3850  }
3851  break;
3852 
3853  case eFourierHalfModeIm:
3854  {
3855  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eFourierSingleModeSpaced);
3856  LibUtilities::BasisKey bkey1(LibUtilities::eFourierHalfModeIm,nummodes_x,pkey1);
3857  returnval.push_back(bkey1);
3858  }
3859  break;
3860 
3861 
3862  case eChebyshev:
3863  {
3864  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eGaussGaussChebyshev);
3865  LibUtilities::BasisKey bkey1(LibUtilities::eChebyshev,nummodes_x,pkey1);
3866  returnval.push_back(bkey1);
3867  }
3868  break;
3869 
3870 
3871 
3872  default:
3873  {
3874  ASSERTL0(false,"Homogeneous expansion can be of Fourier or Chebyshev type only");
3875  }
3876  break;
3877  }
3878 
3879 
3880  switch(type_y)
3881  {
3882  case eFourier:
3883  {
3884  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eFourierEvenlySpaced);
3885  LibUtilities::BasisKey bkey2(LibUtilities::eFourier,nummodes_y,pkey2);
3886  returnval.push_back(bkey2);
3887  }
3888  break;
3889 
3890 
3891  case eFourierSingleMode:
3892  {
3893  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eFourierSingleModeSpaced);
3894  LibUtilities::BasisKey bkey2(LibUtilities::eFourierSingleMode,nummodes_y,pkey2);
3895  returnval.push_back(bkey2);
3896  }
3897  break;
3898 
3899  case eFourierHalfModeRe:
3900  {
3901  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eFourierSingleModeSpaced);
3902  LibUtilities::BasisKey bkey2(LibUtilities::eFourierHalfModeRe,nummodes_y,pkey2);
3903  returnval.push_back(bkey2);
3904  }
3905  break;
3906 
3907  case eFourierHalfModeIm:
3908  {
3909  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eFourierSingleModeSpaced);
3910  LibUtilities::BasisKey bkey2(LibUtilities::eFourierHalfModeIm,nummodes_y,pkey2);
3911  returnval.push_back(bkey2);
3912  }
3913  break;
3914 
3915  case eChebyshev:
3916  {
3917  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eGaussGaussChebyshev);
3918  LibUtilities::BasisKey bkey2(LibUtilities::eChebyshev,nummodes_y,pkey2);
3919  returnval.push_back(bkey2);
3920  }
3921  break;
3922 
3923  default:
3924  {
3925  ASSERTL0(false,"Homogeneous expansion can be of Fourier or Chebyshev type only");
3926  }
3927  break;
3928  }
3929 
3930  switch(type_z)
3931  {
3932  case eFourier:
3933  {
3934  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eFourierEvenlySpaced);
3935  LibUtilities::BasisKey bkey3(LibUtilities::eFourier,nummodes_z,pkey3);
3936  returnval.push_back(bkey3);
3937  }
3938  break;
3939 
3940  case eFourierSingleMode:
3941  {
3942  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eFourierSingleModeSpaced);
3943  LibUtilities::BasisKey bkey3(LibUtilities::eFourierSingleMode,nummodes_z,pkey3);
3944  returnval.push_back(bkey3);
3945  }
3946  break;
3947 
3948  case eFourierHalfModeRe:
3949  {
3950  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eFourierSingleModeSpaced);
3951  LibUtilities::BasisKey bkey3(LibUtilities::eFourierHalfModeRe,nummodes_z,pkey3);
3952  returnval.push_back(bkey3);
3953  }
3954  break;
3955 
3956  case eFourierHalfModeIm:
3957  {
3958  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eFourierSingleModeSpaced);
3959  LibUtilities::BasisKey bkey3(LibUtilities::eFourierHalfModeIm,nummodes_z,pkey3);
3960  returnval.push_back(bkey3);
3961  }
3962  break;
3963 
3964  case eChebyshev:
3965  {
3966  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eGaussGaussChebyshev);
3967  LibUtilities::BasisKey bkey3(LibUtilities::eChebyshev,nummodes_z,pkey3);
3968  returnval.push_back(bkey3);
3969  }
3970  break;
3971 
3972  default:
3973  {
3974  ASSERTL0(false,"Homogeneous expansion can be of Fourier or Chebyshev type only");
3975  }
3976  break;
3977  }
3978  }
3979  break;
3980 
3981  case LibUtilities::eTriangle:
3982  {
3983  ASSERTL0(false,"Homogeneous expansion not defined for this shape");
3984  }
3985  break;
3986 
3987  case LibUtilities::eTetrahedron:
3988  {
3989  ASSERTL0(false,"Homogeneous expansion not defined for this shape");
3990  }
3991  break;
3992 
3993  default:
3994  ASSERTL0(false,"Expansion not defined in switch for this shape");
3995  break;
3996  }
3997 
3998  return returnval;
3999  }
4000 
4001 
4002  /**
4003  *
4004  */
4005  PointGeomSharedPtr MeshGraph::AddVertex(NekDouble x, NekDouble y, NekDouble z)
4006  {
4007  unsigned int nextId = m_vertSet.rbegin()->first + 1;
4008  PointGeomSharedPtr vert(MemoryManager<PointGeom>::AllocateSharedPtr(m_spaceDimension, nextId, x, y, z));
4009  m_vertSet[nextId] = vert;
4010  return vert;
4011  }
4012 
4013  /**
4014  *
4015  */
4016  SegGeomSharedPtr MeshGraph::AddEdge(PointGeomSharedPtr v0, PointGeomSharedPtr v1,
4017  CurveSharedPtr curveDefinition)
4018  {
4019  PointGeomSharedPtr vertices[] = {v0, v1};
4020  SegGeomSharedPtr edge;
4021  int edgeId = m_segGeoms.rbegin()->first + 1;
4022 
4023  if( curveDefinition )
4024  {
4025  edge = MemoryManager<SegGeom>::AllocateSharedPtr(edgeId, m_spaceDimension, vertices, curveDefinition);
4026  }
4027  else
4028  {
4029  edge = MemoryManager<SegGeom>::AllocateSharedPtr(edgeId, m_spaceDimension, vertices);
4030  }
4031  m_segGeoms[edgeId] = edge;
4032  return edge;
4033  }
4034 
4035 
4036  /**
4037  *
4038  */
4039  TriGeomSharedPtr MeshGraph::AddTriangle(SegGeomSharedPtr edges[], StdRegions::Orientation orient[])
4040  {
4041  int indx = m_triGeoms.rbegin()->first + 1;
4042  TriGeomSharedPtr trigeom(MemoryManager<TriGeom>::AllocateSharedPtr(indx, edges, orient));
4043  trigeom->SetGlobalID(indx);
4044 
4045  m_triGeoms[indx] = trigeom;
4046 
4047  return trigeom;
4048  }
4049 
4050 
4051  /**
4052  *
4053  */
4054  QuadGeomSharedPtr MeshGraph::AddQuadrilateral(SegGeomSharedPtr edges[], StdRegions::Orientation orient[])
4055  {
4056  int indx = m_quadGeoms.rbegin()->first + 1;
4057  QuadGeomSharedPtr quadgeom(MemoryManager<QuadGeom>::AllocateSharedPtr(indx, edges, orient));
4058  quadgeom->SetGlobalID(indx);
4059 
4060  m_quadGeoms[indx] = quadgeom;
4061  return quadgeom;
4062  }
4063 
4064 
4065  /**
4066  *
4067  */
4068  PrismGeomSharedPtr MeshGraph::AddPrism(TriGeomSharedPtr tfaces[PrismGeom::kNtfaces],
4069  QuadGeomSharedPtr qfaces[PrismGeom::kNqfaces])
4070  {
4071  // Setting the orientation is disabled in the reader. Why?
4072  Geometry2DSharedPtr faces[] = { qfaces[0], tfaces[0], qfaces[1], tfaces[1], qfaces[2] };
4073  unsigned int index = m_prismGeoms.rbegin()->first + 1;
4074  PrismGeomSharedPtr prismgeom(MemoryManager<PrismGeom>::AllocateSharedPtr(faces));
4075  prismgeom->SetGlobalID(index);
4076 
4077  m_prismGeoms[index] = prismgeom;
4078  return prismgeom;
4079  }
4080 
4081 
4082  /**
4083  *
4084  */
4085  TetGeomSharedPtr MeshGraph::AddTetrahedron(TriGeomSharedPtr tfaces[TetGeom::kNtfaces])
4086  {
4087  unsigned int index = m_tetGeoms.rbegin()->first + 1;
4088  TetGeomSharedPtr tetgeom(MemoryManager<TetGeom>::AllocateSharedPtr(tfaces));
4089  tetgeom->SetGlobalID(index);
4090 
4091  m_tetGeoms[index] = tetgeom;
4092  return tetgeom;
4093  }
4094 
4095 
4096  /**
4097  *
4098  */
4099  PyrGeomSharedPtr MeshGraph::AddPyramid(TriGeomSharedPtr tfaces[PyrGeom::kNtfaces],
4100  QuadGeomSharedPtr qfaces[PyrGeom::kNqfaces])
4101  {
4102  Geometry2DSharedPtr faces[] = { qfaces[0], tfaces[0], tfaces[1], tfaces[2], tfaces[3] };
4103  unsigned int index = m_pyrGeoms.rbegin()->first + 1;
4104 
4105  PyrGeomSharedPtr pyrgeom(MemoryManager<PyrGeom>::AllocateSharedPtr(faces));
4106  pyrgeom->SetGlobalID(index);
4107 
4108  m_pyrGeoms[index] = pyrgeom;
4109  return pyrgeom;
4110  }
4111 
4112 
4113  /**
4114  *
4115  */
4116  HexGeomSharedPtr MeshGraph::AddHexahedron(QuadGeomSharedPtr qfaces[HexGeom::kNqfaces])
4117  {
4118  unsigned int index = m_hexGeoms.rbegin()->first + 1;
4119  HexGeomSharedPtr hexgeom(MemoryManager<HexGeom>::AllocateSharedPtr(qfaces));
4120  hexgeom->SetGlobalID(index);
4121  m_hexGeoms[index] = hexgeom;
4122  return hexgeom;
4123  }
4124 
4125 
4126  /**
4127  * Generate a single vector of Expansion structs mapping global element
4128  * ID to a corresponding Geometry shared pointer and basis key.
4129  *
4130  * Expansion map ensures elements which appear in multiple composites
4131  * within the domain are only listed once.
4132  */
4133  ExpansionMapShPtr MeshGraph::SetUpExpansionMap(void)
4134  {
4135  ExpansionMapShPtr returnval;
4136  returnval = MemoryManager<ExpansionMap>::AllocateSharedPtr();
4137 
4138  for(int d = 0; d < m_domain.size(); ++d)
4139  {
4140  CompositeMap::const_iterator compIter;
4141 
4142  for (compIter = m_domain[d].begin(); compIter != m_domain[d].end(); ++compIter)
4143  {
4144  GeometryVector::const_iterator x;
4145  for (x = compIter->second->begin(); x != compIter->second->end(); ++x)
4146  {
4147  LibUtilities::BasisKeyVector def;
4148  ExpansionShPtr expansionElementShPtr =
4149  MemoryManager<Expansion>::AllocateSharedPtr(*x, def);
4150  int id = (*x)->GetGlobalID();
4151  (*returnval)[id] = expansionElementShPtr;
4152  }
4153  }
4154  }
4155 
4156  return returnval;
4157  }
4158  }; //end of namespace
4159 }; //end of namespace
StdPrismExp.h
Nektar::SpatialDomains::MeshGraph::SetExpansions
void SetExpansions(std::vector< LibUtilities::FieldDefinitionsSharedPtr > &fielddef)
Sets expansions given field definitions.
Definition: MeshGraph.cpp:2363
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:3105
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:4068
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:4005
Nektar::SpatialDomains::MeshGraph::AddQuadrilateral
QuadGeomSharedPtr AddQuadrilateral(SegGeomSharedPtr edges[], StdRegions::Orientation orient[])
Definition: MeshGraph.cpp:4054
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:3206
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:3138
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:4039
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:1971
Nektar::SpatialDomains::MeshGraph::GetCompositeItem
GeometrySharedPtr GetCompositeItem(int whichComposite, int whichItem)
Definition: MeshGraph.cpp:2233
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:4016
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:4099
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:3798
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:2271
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:1985
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:2023
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:3185
Nektar::SpatialDomains::MeshGraph::GetExpansion
ExpansionShPtr GetExpansion(GeometrySharedPtr geom, const std::string variable="DefaultVar")
Definition: MeshGraph.cpp:2346
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:4116
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:4133
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:3063
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:4085
StdTetExp.h
Nektar::SpatialDomains::ExpansionMapIter
std::map< int, ExpansionShPtr >::iterator ExpansionMapIter
Definition: MeshGraph.h:175
Generated on Wed Mar 8 2017 12:28:57 for Nektar++ by   doxygen 1.8.8