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