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