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