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  ExpansionMapShPtr expansionMap = m_expansionMapShPtrMap.find(variable)->second;
2329 
2330  iter = expansionMap->find(geom->GetGlobalID());
2331  ASSERTL1(iter != expansionMap->end(),
2332  "Could not find expansion " +
2333  boost::lexical_cast<string>(geom->GetGlobalID()) +
2334  " in expansion for variable " + variable);
2335  return iter->second;
2336  }
2337 
2338 
2339  /**
2340  *
2341  */
2342  void MeshGraph::SetExpansions(
2343  std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
2344  {
2345  int i, j, k, cnt, id;
2346  GeometrySharedPtr geom;
2347 
2348  ExpansionMapShPtr expansionMap;
2349 
2350  // Loop over fields and determine unique fields string and
2351  // declare whole expansion list
2352  for(i = 0; i < fielddef.size(); ++i)
2353  {
2354  for(j = 0; j < fielddef[i]->m_fields.size(); ++j)
2355  {
2356  std::string field = fielddef[i]->m_fields[j];
2357  if(m_expansionMapShPtrMap.count(field) == 0)
2358  {
2359  expansionMap = MemoryManager<ExpansionMap>::AllocateSharedPtr();
2360  m_expansionMapShPtrMap[field] = expansionMap;
2361 
2362  // check to see if DefaultVar also not set and
2363  // if so assign it to this expansion
2364  if(m_expansionMapShPtrMap.count("DefaultVar") == 0)
2365  {
2366  m_expansionMapShPtrMap["DefaultVar"] = expansionMap;
2367  }
2368 
2369  // loop over all elements and set expansion
2370  for(k = 0; k < fielddef.size(); ++k)
2371  {
2372  for(int h = 0; h < fielddef[k]->m_fields.size(); ++h)
2373  {
2374  if(fielddef[k]->m_fields[h] == field)
2375  {
2376  expansionMap = m_expansionMapShPtrMap.find(field)->second;
2377  LibUtilities::BasisKeyVector def;
2378 
2379  for(int g = 0; g < fielddef[k]->m_elementIDs.size(); ++g)
2380  {
2381  ExpansionShPtr tmpexp =
2382  MemoryManager<Expansion>::AllocateSharedPtr(geom, def);
2383  (*expansionMap)[fielddef[k]->m_elementIDs[g]] = tmpexp;
2384  }
2385  }
2386  }
2387  }
2388  }
2389  }
2390  }
2391 
2392  // loop over all elements find the geometry shared ptr and
2393  // set up basiskey vector
2394  for(i = 0; i < fielddef.size(); ++i)
2395  {
2396  cnt = 0;
2397  std::vector<std::string> fields = fielddef[i]->m_fields;
2398  std::vector<unsigned int> nmodes = fielddef[i]->m_numModes;
2399  std::vector<LibUtilities::BasisType> basis = fielddef[i]->m_basis;
2400  bool pointDef = fielddef[i]->m_pointsDef;
2401  bool numPointDef = fielddef[i]->m_numPointsDef;
2402 
2403  // Check points and numpoints
2404  std::vector<unsigned int> npoints = fielddef[i]->m_numPoints;
2405  std::vector<LibUtilities::PointsType> points = fielddef[i]->m_points;
2406 
2407  bool UniOrder = fielddef[i]->m_uniOrder;
2408 
2409  for (j = 0; j < fielddef[i]->m_elementIDs.size(); ++j)
2410  {
2411 
2412  LibUtilities::BasisKeyVector bkeyvec;
2413  id = fielddef[i]->m_elementIDs[j];
2414 
2415  switch (fielddef[i]->m_shapeType)
2416  {
2417  case LibUtilities::eSegment:
2418  {
2419  if(m_segGeoms.count(fielddef[i]->m_elementIDs[j]) == 0)
2420  {
2421  // skip element likely from parallel read
2422  continue;
2423  }
2424  geom = m_segGeoms[fielddef[i]->m_elementIDs[j]];
2425 
2426  LibUtilities::PointsKey pkey(nmodes[cnt]+1, LibUtilities::eGaussLobattoLegendre);
2427 
2428  if(numPointDef&&pointDef)
2429  {
2430  const LibUtilities::PointsKey pkey1(npoints[cnt], points[0]);
2431  pkey = pkey1;
2432  }
2433  else if(!numPointDef&&pointDef)
2434  {
2435  const LibUtilities::PointsKey pkey1(nmodes[cnt]+1, points[0]);
2436  pkey = pkey1;
2437  }
2438  else if(numPointDef&&!pointDef)
2439  {
2440  const LibUtilities::PointsKey pkey1(npoints[cnt], LibUtilities::eGaussLobattoLegendre);
2441  pkey = pkey1;
2442  }
2443 
2444  LibUtilities::BasisKey bkey(basis[0], nmodes[cnt], pkey);
2445 
2446  if(!UniOrder)
2447  {
2448  cnt++;
2449  cnt += fielddef[i]->m_numHomogeneousDir;
2450  }
2451  bkeyvec.push_back(bkey);
2452  }
2453  break;
2454  case LibUtilities::eTriangle:
2455  {
2456  if(m_triGeoms.count(fielddef[i]->m_elementIDs[j]) == 0)
2457  {
2458  // skip element likely from parallel read
2459  continue;
2460  }
2461  geom = m_triGeoms[fielddef[i]->m_elementIDs[j]];
2462 
2463  LibUtilities::PointsKey pkey(nmodes[cnt]+1, LibUtilities::eGaussLobattoLegendre);
2464  if(numPointDef&&pointDef)
2465  {
2466  const LibUtilities::PointsKey pkey2(npoints[cnt], points[0]);
2467  pkey = pkey2;
2468  }
2469  else if(!numPointDef&&pointDef)
2470  {
2471  const LibUtilities::PointsKey pkey2(nmodes[cnt]+1, points[0]);
2472  pkey = pkey2;
2473  }
2474  else if(numPointDef&&!pointDef)
2475  {
2476  const LibUtilities::PointsKey pkey2(npoints[cnt], LibUtilities::eGaussLobattoLegendre);
2477  pkey = pkey2;
2478  }
2479  LibUtilities::BasisKey bkey(basis[0], nmodes[cnt], pkey);
2480 
2481  bkeyvec.push_back(bkey);
2482 
2483  LibUtilities::PointsKey pkey1(nmodes[cnt+1], LibUtilities::eGaussRadauMAlpha1Beta0);
2484  if(numPointDef&&pointDef)
2485  {
2486  const LibUtilities::PointsKey pkey2(npoints[cnt+1], points[1]);
2487  pkey1 = pkey2;
2488  }
2489  else if(!numPointDef&&pointDef)
2490  {
2491  const LibUtilities::PointsKey pkey2(nmodes[cnt+1], points[1]);
2492  pkey1 = pkey2;
2493  }
2494  else if(numPointDef&&!pointDef)
2495  {
2496  const LibUtilities::PointsKey pkey2(npoints[cnt+1], LibUtilities::eGaussRadauMAlpha1Beta0);
2497  pkey1 = pkey2;
2498  }
2499  LibUtilities::BasisKey bkey1(basis[1], nmodes[cnt+1], pkey1);
2500  bkeyvec.push_back(bkey1);
2501 
2502  if(!UniOrder)
2503  {
2504  cnt += 2;
2505  cnt += fielddef[i]->m_numHomogeneousDir;
2506  }
2507  }
2508  break;
2509  case LibUtilities::eQuadrilateral:
2510  {
2511  if(m_quadGeoms.count(fielddef[i]->m_elementIDs[j]) == 0)
2512  {
2513  // skip element likely from parallel read
2514  continue;
2515  }
2516 
2517  geom = m_quadGeoms[fielddef[i]->m_elementIDs[j]];
2518 
2519  for(int b = 0; b < 2; ++b)
2520  {
2521  LibUtilities::PointsKey pkey(nmodes[cnt+b]+1, LibUtilities::eGaussLobattoLegendre);
2522 
2523  if(numPointDef&&pointDef)
2524  {
2525  const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
2526  pkey = pkey2;
2527  }
2528  else if(!numPointDef&&pointDef)
2529  {
2530  const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
2531  pkey = pkey2;
2532  }
2533  else if(numPointDef&&!pointDef)
2534  {
2535  const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
2536  pkey = pkey2;
2537  }
2538  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2539  bkeyvec.push_back(bkey);
2540  }
2541 
2542  if(!UniOrder)
2543  {
2544  cnt += 2;
2545  cnt += fielddef[i]->m_numHomogeneousDir;
2546  }
2547  }
2548  break;
2549 
2550  case LibUtilities::eTetrahedron:
2551  {
2552  k = fielddef[i]->m_elementIDs[j];
2553 
2554  // allow for possibility that fielddef is
2555  // larger than m_graph which can happen in
2556  // parallel runs
2557  if(m_tetGeoms.count(k) == 0)
2558  {
2559  continue;
2560  }
2561  geom = m_tetGeoms[k];
2562 
2563  {
2564  LibUtilities::PointsKey pkey(nmodes[cnt]+1, LibUtilities::eGaussLobattoLegendre);
2565 
2566  if(numPointDef&&pointDef)
2567  {
2568  const LibUtilities::PointsKey pkey2(npoints[cnt],points[0]);
2569  pkey = pkey2;
2570  }
2571  else if(!numPointDef&&pointDef)
2572  {
2573  const LibUtilities::PointsKey pkey2(nmodes[cnt]+1,points[0]);
2574  pkey = pkey2;
2575  }
2576  else if(numPointDef&&!pointDef)
2577  {
2578  const LibUtilities::PointsKey pkey2(npoints[cnt],LibUtilities::eGaussLobattoLegendre);
2579  pkey = pkey2;
2580  }
2581 
2582  LibUtilities::BasisKey bkey(basis[0],nmodes[cnt],pkey);
2583 
2584  bkeyvec.push_back(bkey);
2585  }
2586  {
2587  LibUtilities::PointsKey pkey(nmodes[cnt+1], LibUtilities::eGaussRadauMAlpha1Beta0);
2588 
2589  if(numPointDef&&pointDef)
2590  {
2591  const LibUtilities::PointsKey pkey2(npoints[cnt+1],points[1]);
2592  pkey = pkey2;
2593  }
2594  else if(!numPointDef&&pointDef)
2595  {
2596  const LibUtilities::PointsKey pkey2(nmodes[cnt+1]+1,points[1]);
2597  pkey = pkey2;
2598  }
2599  else if(numPointDef&&!pointDef)
2600  {
2601  const LibUtilities::PointsKey pkey2(npoints[cnt+1],LibUtilities::eGaussRadauMAlpha1Beta0);
2602  pkey = pkey2;
2603  }
2604 
2605  LibUtilities::BasisKey bkey(basis[1],nmodes[cnt+1],pkey);
2606 
2607  bkeyvec.push_back(bkey);
2608  }
2609 
2610  {
2611  LibUtilities::PointsKey pkey(nmodes[cnt+2], LibUtilities::eGaussRadauMAlpha2Beta0);
2612 
2613  if(numPointDef&&pointDef)
2614  {
2615  const LibUtilities::PointsKey pkey2(npoints[cnt+2],points[2]);
2616  pkey = pkey2;
2617  }
2618  else if(!numPointDef&&pointDef)
2619  {
2620  const LibUtilities::PointsKey pkey2(nmodes[cnt+2]+1,points[2]);
2621  pkey = pkey2;
2622  }
2623  else if(numPointDef&&!pointDef)
2624  {
2625  const LibUtilities::PointsKey pkey2(npoints[cnt+2],LibUtilities::eGaussRadauMAlpha1Beta0);
2626  pkey = pkey2;
2627  }
2628 
2629  LibUtilities::BasisKey bkey(basis[2],nmodes[cnt+2],pkey);
2630 
2631  bkeyvec.push_back(bkey);
2632  }
2633 
2634  if(!UniOrder)
2635  {
2636  cnt += 3;
2637  }
2638  }
2639  break;
2640  case LibUtilities::ePrism:
2641  {
2642  k = fielddef[i]->m_elementIDs[j];
2643  if(m_prismGeoms.count(k) == 0)
2644  {
2645  continue;
2646  }
2647  geom = m_prismGeoms[k];
2648 
2649  for(int b = 0; b < 2; ++b)
2650  {
2651  LibUtilities::PointsKey pkey(nmodes[cnt+b]+1,LibUtilities::eGaussLobattoLegendre);
2652 
2653  if(numPointDef&&pointDef)
2654  {
2655  const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
2656  pkey = pkey2;
2657  }
2658  else if(!numPointDef&&pointDef)
2659  {
2660  const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
2661  pkey = pkey2;
2662  }
2663  else if(numPointDef&&!pointDef)
2664  {
2665  const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
2666  pkey = pkey2;
2667  }
2668 
2669  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2670  bkeyvec.push_back(bkey);
2671  }
2672 
2673  {
2674  LibUtilities::PointsKey pkey(nmodes[cnt+2],LibUtilities::eGaussRadauMAlpha1Beta0);
2675 
2676  if(numPointDef&&pointDef)
2677  {
2678  const LibUtilities::PointsKey pkey2(npoints[cnt+2],points[2]);
2679  pkey = pkey2;
2680  }
2681  else if(!numPointDef&&pointDef)
2682  {
2683  const LibUtilities::PointsKey pkey2(nmodes[cnt+2]+1,points[2]);
2684  pkey = pkey2;
2685  }
2686  else if(numPointDef&&!pointDef)
2687  {
2688  const LibUtilities::PointsKey pkey2(npoints[cnt+2],LibUtilities::eGaussLobattoLegendre);
2689  pkey = pkey2;
2690  }
2691 
2692  LibUtilities::BasisKey bkey(basis[2],nmodes[cnt+2],pkey);
2693  bkeyvec.push_back(bkey);
2694  }
2695 
2696  if(!UniOrder)
2697  {
2698  cnt += 3;
2699  }
2700  }
2701  break;
2702  case LibUtilities::ePyramid:
2703  {
2704  k = fielddef[i]->m_elementIDs[j];
2705  ASSERTL0(m_pyrGeoms.find(k) != m_pyrGeoms.end(),
2706  "Failed to find geometry with same global id");
2707  geom = m_pyrGeoms[k];
2708 
2709 
2710  for(int b = 0; b < 2; ++b)
2711  {
2712  LibUtilities::PointsKey pkey(nmodes[cnt+b]+1,LibUtilities::eGaussLobattoLegendre);
2713 
2714  if(numPointDef&&pointDef)
2715  {
2716  const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
2717  pkey = pkey2;
2718  }
2719  else if(!numPointDef&&pointDef)
2720  {
2721  const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
2722  pkey = pkey2;
2723  }
2724  else if(numPointDef&&!pointDef)
2725  {
2726  const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
2727  pkey = pkey2;
2728  }
2729 
2730  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2731  bkeyvec.push_back(bkey);
2732  }
2733 
2734  {
2735  LibUtilities::PointsKey pkey(nmodes[cnt+2],LibUtilities::eGaussRadauMAlpha2Beta0);
2736 
2737  if(numPointDef&&pointDef)
2738  {
2739  const LibUtilities::PointsKey pkey2(npoints[cnt+2],points[2]);
2740  pkey = pkey2;
2741  }
2742  else if(!numPointDef&&pointDef)
2743  {
2744  const LibUtilities::PointsKey pkey2(nmodes[cnt+2]+1,points[2]);
2745  pkey = pkey2;
2746  }
2747  else if(numPointDef&&!pointDef)
2748  {
2749  const LibUtilities::PointsKey pkey2(npoints[cnt+2],LibUtilities::eGaussLobattoLegendre);
2750  pkey = pkey2;
2751  }
2752 
2753  LibUtilities::BasisKey bkey(basis[2],nmodes[cnt+2],pkey);
2754  bkeyvec.push_back(bkey);
2755  }
2756 
2757  if(!UniOrder)
2758  {
2759  cnt += 3;
2760  }
2761  }
2762  break;
2763  case LibUtilities::eHexahedron:
2764  {
2765  k = fielddef[i]->m_elementIDs[j];
2766  if(m_hexGeoms.count(k) == 0)
2767  {
2768  continue;
2769  }
2770 
2771  geom = m_hexGeoms[k];
2772 
2773  for(int b = 0; b < 3; ++b)
2774  {
2775  LibUtilities::PointsKey pkey(nmodes[cnt+b],LibUtilities::eGaussLobattoLegendre);
2776 
2777  if(numPointDef&&pointDef)
2778  {
2779  const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
2780  pkey = pkey2;
2781  }
2782  else if(!numPointDef&&pointDef)
2783  {
2784  const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
2785  pkey = pkey2;
2786  }
2787  else if(numPointDef&&!pointDef)
2788  {
2789  const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
2790  pkey = pkey2;
2791  }
2792 
2793  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2794  bkeyvec.push_back(bkey);
2795  }
2796 
2797  if(!UniOrder)
2798  {
2799  cnt += 3;
2800  }
2801  }
2802  break;
2803  default:
2804  ASSERTL0(false,"Need to set up for pyramid and prism 3D Expansions");
2805  break;
2806  }
2807 
2808  for(k = 0; k < fields.size(); ++k)
2809  {
2810  expansionMap = m_expansionMapShPtrMap.find(fields[k])->second;
2811  if((*expansionMap).find(id) != (*expansionMap).end())
2812  {
2813  (*expansionMap)[id]->m_geomShPtr = geom;
2814  (*expansionMap)[id]->m_basisKeyVector = bkeyvec;
2815  }
2816  }
2817  }
2818  }
2819  }
2820 
2821 
2822  /**
2823  *
2824  */
2825  void MeshGraph::SetExpansions(
2826  std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef,
2827  std::vector< std::vector<LibUtilities::PointsType> > &pointstype)
2828  {
2829  int i,j,k,g,h,cnt,id;
2830  GeometrySharedPtr geom;
2831 
2832  ExpansionMapShPtr expansionMap;
2833 
2834  // Loop over fields and determine unique fields string and
2835  // declare whole expansion list
2836  for(i = 0; i < fielddef.size(); ++i)
2837  {
2838  for(j = 0; j < fielddef[i]->m_fields.size(); ++j)
2839  {
2840  std::string field = fielddef[i]->m_fields[j];
2841  if(m_expansionMapShPtrMap.count(field) == 0)
2842  {
2843  expansionMap = MemoryManager<ExpansionMap>::AllocateSharedPtr();
2844  m_expansionMapShPtrMap[field] = expansionMap;
2845 
2846  // check to see if DefaultVar also not set and
2847  // if so assign it to this expansion
2848  if(m_expansionMapShPtrMap.count("DefaultVar") == 0)
2849  {
2850  m_expansionMapShPtrMap["DefaultVar"] = expansionMap;
2851  }
2852 
2853  // loop over all elements and set expansion
2854  for(k = 0; k < fielddef.size(); ++k)
2855  {
2856  for(h = 0; h < fielddef[k]->m_fields.size(); ++h)
2857  {
2858  if(fielddef[k]->m_fields[h] == field)
2859  {
2860  expansionMap = m_expansionMapShPtrMap.find(field)->second;
2861  LibUtilities::BasisKeyVector def;
2862 
2863  for(g = 0; g < fielddef[k]->m_elementIDs.size(); ++g)
2864  {
2865  ExpansionShPtr tmpexp =
2866  MemoryManager<Expansion>::AllocateSharedPtr(geom, def);
2867  (*expansionMap)[fielddef[k]->m_elementIDs[g]] = tmpexp;
2868  }
2869  }
2870  }
2871  }
2872  }
2873  }
2874  }
2875 
2876 
2877  // loop over all elements find the geometry shared ptr and
2878  // set up basiskey vector
2879  for(i = 0; i < fielddef.size(); ++i)
2880  {
2881  cnt = 0;
2882  std::vector<std::string> fields = fielddef[i]->m_fields;
2883  std::vector<unsigned int> nmodes = fielddef[i]->m_numModes;
2884  std::vector<LibUtilities::BasisType> basis = fielddef[i]->m_basis;
2885  bool UniOrder = fielddef[i]->m_uniOrder;
2886 
2887  for(j = 0; j < fielddef[i]->m_elementIDs.size(); ++j)
2888  {
2889  LibUtilities::BasisKeyVector bkeyvec;
2890  id = fielddef[i]->m_elementIDs[j];
2891 
2892  switch(fielddef[i]->m_shapeType)
2893  {
2894  case LibUtilities::eSegment:
2895  {
2896  k = fielddef[i]->m_elementIDs[j];
2897  ASSERTL0(m_segGeoms.find(k) != m_segGeoms.end(),
2898  "Failed to find geometry with same global id.");
2899  geom = m_segGeoms[k];
2900 
2901  const LibUtilities::PointsKey pkey(nmodes[cnt], pointstype[i][0]);
2902  LibUtilities::BasisKey bkey(basis[0], nmodes[cnt], pkey);
2903  if(!UniOrder)
2904  {
2905  cnt++;
2906  }
2907  bkeyvec.push_back(bkey);
2908  }
2909  break;
2910  case LibUtilities::eTriangle:
2911  {
2912  k = fielddef[i]->m_elementIDs[j];
2913  ASSERTL0(m_triGeoms.find(k) != m_triGeoms.end(),
2914  "Failed to find geometry with same global id.");
2915  geom = m_triGeoms[k];
2916  for(int b = 0; b < 2; ++b)
2917  {
2918  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
2919  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2920  bkeyvec.push_back(bkey);
2921  }
2922 
2923  if(!UniOrder)
2924  {
2925  cnt += 2;
2926  }
2927  }
2928  break;
2929  case LibUtilities::eQuadrilateral:
2930  {
2931  k = fielddef[i]->m_elementIDs[j];
2932  ASSERTL0(m_quadGeoms.find(k) != m_quadGeoms.end(),
2933  "Failed to find geometry with same global id");
2934  geom = m_quadGeoms[k];
2935 
2936  for(int b = 0; b < 2; ++b)
2937  {
2938  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
2939  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2940  bkeyvec.push_back(bkey);
2941  }
2942 
2943  if(!UniOrder)
2944  {
2945  cnt += 2;
2946  }
2947  }
2948  break;
2949  case LibUtilities::eTetrahedron:
2950  {
2951  k = fielddef[i]->m_elementIDs[j];
2952  ASSERTL0(m_tetGeoms.find(k) != m_tetGeoms.end(),
2953  "Failed to find geometry with same global id");
2954  geom = m_tetGeoms[k];
2955 
2956  for(int b = 0; b < 3; ++b)
2957  {
2958  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
2959  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2960  bkeyvec.push_back(bkey);
2961  }
2962 
2963  if(!UniOrder)
2964  {
2965  cnt += 2;
2966  }
2967  }
2968  break;
2969  case LibUtilities::ePyramid:
2970  {
2971  k = fielddef[i]->m_elementIDs[j];
2972  ASSERTL0(m_pyrGeoms.find(k) != m_pyrGeoms.end(),
2973  "Failed to find geometry with same global id");
2974  geom = m_pyrGeoms[k];
2975 
2976  for(int b = 0; b < 3; ++b)
2977  {
2978  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
2979  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2980  bkeyvec.push_back(bkey);
2981  }
2982 
2983  if(!UniOrder)
2984  {
2985  cnt += 2;
2986  }
2987  }
2988  break;
2989  case LibUtilities::ePrism:
2990  {
2991  k = fielddef[i]->m_elementIDs[j];
2992  ASSERTL0(m_prismGeoms.find(k) != m_prismGeoms.end(),
2993  "Failed to find geometry with same global id");
2994  geom = m_prismGeoms[k];
2995 
2996  for(int b = 0; b < 3; ++b)
2997  {
2998  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
2999  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
3000  bkeyvec.push_back(bkey);
3001  }
3002 
3003  if(!UniOrder)
3004  {
3005  cnt += 2;
3006  }
3007  }
3008  break;
3009  case LibUtilities::eHexahedron:
3010  {
3011  k = fielddef[i]->m_elementIDs[j];
3012  ASSERTL0(m_hexGeoms.find(k) != m_hexGeoms.end(),
3013  "Failed to find geometry with same global id");
3014  geom = m_hexGeoms[k];
3015 
3016  for(int b = 0; b < 3; ++b)
3017  {
3018  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
3019  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
3020  bkeyvec.push_back(bkey);
3021  }
3022 
3023  if(!UniOrder)
3024  {
3025  cnt += 2;
3026  }
3027  }
3028  break;
3029  default:
3030  ASSERTL0(false,"Need to set up for pyramid and prism 3D Expansions");
3031  break;
3032  }
3033 
3034  for(k = 0; k < fields.size(); ++k)
3035  {
3036  expansionMap = m_expansionMapShPtrMap.find(fields[k])->second;
3037  if((*expansionMap).find(id) != (*expansionMap).end())
3038  {
3039  (*expansionMap)[id]->m_geomShPtr = geom;
3040  (*expansionMap)[id]->m_basisKeyVector = bkeyvec;
3041  }
3042  }
3043  }
3044  }
3045  }
3046 
3047  /**
3048  * \brief Reset all points keys to have equispaced points with
3049  * optional arguemt of \a npoints which redefines how many
3050  * points are to be used.
3051  */
3052  void MeshGraph::SetExpansionsToEvenlySpacedPoints(int npoints)
3053  {
3054  ExpansionMapShPtrMapIter it;
3055 
3056  // iterate over all defined expansions
3057  for(it = m_expansionMapShPtrMap.begin(); it != m_expansionMapShPtrMap.end(); ++it)
3058  {
3059  ExpansionMapIter expIt;
3060 
3061  for(expIt = it->second->begin(); expIt != it->second->end(); ++expIt)
3062  {
3063  for(int i = 0; i < expIt->second->m_basisKeyVector.size(); ++i)
3064  {
3065  LibUtilities::BasisKey bkeyold = expIt->second->m_basisKeyVector[i];
3066 
3067  int npts;
3068 
3069  if(npoints) // use input
3070  {
3071  npts = npoints;
3072  }
3073  else
3074  {
3075  npts = bkeyold.GetNumModes();
3076  }
3077 
3078 
3079  const LibUtilities::PointsKey pkey(npts,LibUtilities::ePolyEvenlySpaced);
3080  LibUtilities::BasisKey bkeynew(bkeyold.GetBasisType(),bkeyold.GetNumModes(), pkey);
3081  expIt->second->m_basisKeyVector[i] = bkeynew;
3082 
3083  }
3084  }
3085  }
3086  }
3087 
3088  /**
3089  * \brief Reset all points keys to have expansion order of \a
3090  * nmodes. we keep the point distribution the same and make
3091  * the number of points the same difference from the number
3092  * of modes as the original expansion definition
3093  */
3094  void MeshGraph::SetExpansionsToPolyOrder(int nmodes)
3095  {
3096  ExpansionMapShPtrMapIter it;
3097 
3098  // iterate over all defined expansions
3099  for(it = m_expansionMapShPtrMap.begin(); it != m_expansionMapShPtrMap.end(); ++it)
3100  {
3101  ExpansionMapIter expIt;
3102 
3103  for(expIt = it->second->begin(); expIt != it->second->end(); ++expIt)
3104  {
3105  for(int i = 0; i < expIt->second->m_basisKeyVector.size(); ++i)
3106  {
3107  LibUtilities::BasisKey bkeyold = expIt->second->m_basisKeyVector[i];
3108 
3109  int npts = nmodes + (bkeyold.GetNumPoints() - bkeyold.GetNumModes());
3110 
3111  const LibUtilities::PointsKey pkey(npts,bkeyold.GetPointsType());
3112  LibUtilities::BasisKey bkeynew(bkeyold.GetBasisType(),nmodes, pkey);
3113  expIt->second->m_basisKeyVector[i] = bkeynew;
3114 
3115  }
3116  }
3117  }
3118  }
3119 
3120 
3121  /**
3122  * \brief Reset all points keys to have expansion order of \a
3123  * nmodes. we keep the point distribution the same and make
3124  * the number of points the same difference from the number
3125  * of modes as the original expansion definition
3126  */
3127  void MeshGraph::SetExpansionsToPointOrder(int npts)
3128  {
3129  ExpansionMapShPtrMapIter it;
3130 
3131  // iterate over all defined expansions
3132  for (it = m_expansionMapShPtrMap.begin();
3133  it != m_expansionMapShPtrMap.end();
3134  ++it)
3135  {
3136  ExpansionMapIter expIt;
3137 
3138  for (expIt = it->second->begin();
3139  expIt != it->second->end();
3140  ++expIt)
3141  {
3142  for(int i = 0;
3143  i < expIt->second->m_basisKeyVector.size();
3144  ++i)
3145  {
3146  LibUtilities::BasisKey bkeyold =
3147  expIt->second->m_basisKeyVector[i];
3148 
3149  const LibUtilities::PointsKey pkey(
3150  npts, bkeyold.GetPointsType());
3151 
3152  LibUtilities::BasisKey bkeynew(bkeyold.GetBasisType(),
3153  bkeyold.GetNumModes(),
3154  pkey);
3155  expIt->second->m_basisKeyVector[i] = bkeynew;
3156  }
3157  }
3158  }
3159  }
3160 
3161 
3162  /**
3163  * For each element of shape given by \a shape in field \a
3164  * var, replace the current BasisKeyVector describing the
3165  * expansion in each dimension, with the one provided by \a
3166  * keys.
3167  *
3168  * @TODO: Allow selection of elements through a CompositeVector,
3169  * as well as by type.
3170  *
3171  * @param shape The shape of elements to be changed.
3172  * @param keys The new basis vector to apply to those elements.
3173  */
3174  void MeshGraph::SetBasisKey(LibUtilities::ShapeType shape,
3175  LibUtilities::BasisKeyVector &keys,
3176  std::string var)
3177  {
3178  ExpansionMapIter elemIter;
3179 
3180  ExpansionMapShPtr expansionMap = m_expansionMapShPtrMap.find(var)->second;
3181 
3182  for (elemIter = expansionMap->begin(); elemIter != expansionMap->end(); ++elemIter)
3183  {
3184  if ((elemIter->second)->m_geomShPtr->GetShapeType() == shape)
3185  {
3186  (elemIter->second)->m_basisKeyVector = keys;
3187  }
3188  }
3189  }
3190 
3191 
3192  /**
3193  *
3194  */
3195  LibUtilities::BasisKeyVector MeshGraph::DefineBasisKeyFromExpansionType(
3196  GeometrySharedPtr in,
3197  ExpansionType type,
3198  const int nummodes)
3199  {
3200  LibUtilities::BasisKeyVector returnval;
3201 
3202  LibUtilities::ShapeType shape= in->GetShapeType();
3203 
3204  int quadoffset = 1;
3205  switch(type)
3206  {
3207  case eModified:
3208  case eModifiedGLLRadau10:
3209  quadoffset = 1;
3210  break;
3211  case eModifiedQuadPlus1:
3212  quadoffset = 2;
3213  break;
3214  case eModifiedQuadPlus2:
3215  quadoffset = 3;
3216  break;
3217  default:
3218  break;
3219  }
3220 
3221  switch(type)
3222  {
3223  case eModified:
3224  case eModifiedQuadPlus1:
3225  case eModifiedQuadPlus2:
3226  case eModifiedGLLRadau10:
3227  {
3228  switch (shape)
3229  {
3230  case LibUtilities::eSegment:
3231  {
3232  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3233  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3234  returnval.push_back(bkey);
3235  }
3236  break;
3237  case LibUtilities::eQuadrilateral:
3238  {
3239  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3240  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3241  returnval.push_back(bkey);
3242  returnval.push_back(bkey);
3243  }
3244  break;
3245  case LibUtilities::eHexahedron:
3246  {
3247  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3248  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3249  returnval.push_back(bkey);
3250  returnval.push_back(bkey);
3251  returnval.push_back(bkey);
3252  }
3253  break;
3254  case LibUtilities::eTriangle:
3255  {
3256  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3257  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3258  returnval.push_back(bkey);
3259 
3260  const LibUtilities::PointsKey pkey1(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha1Beta0);
3261  LibUtilities::BasisKey bkey1(LibUtilities::eModified_B, nummodes, pkey1);
3262 
3263  returnval.push_back(bkey1);
3264  }
3265  break;
3266  case LibUtilities::eTetrahedron:
3267  {
3268  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3269  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3270  returnval.push_back(bkey);
3271 
3272  const LibUtilities::PointsKey pkey1(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha1Beta0);
3273  LibUtilities::BasisKey bkey1(LibUtilities::eModified_B, nummodes, pkey1);
3274  returnval.push_back(bkey1);
3275 
3276  if(type == eModifiedGLLRadau10)
3277  {
3278  const LibUtilities::PointsKey pkey2(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha1Beta0);
3279  LibUtilities::BasisKey bkey2(LibUtilities::eModified_C, nummodes, pkey2);
3280  returnval.push_back(bkey2);
3281  }
3282  else
3283  {
3284  const LibUtilities::PointsKey pkey2(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha2Beta0);
3285  LibUtilities::BasisKey bkey2(LibUtilities::eModified_C, nummodes, pkey2);
3286  returnval.push_back(bkey2);
3287  }
3288  }
3289  break;
3290  case LibUtilities::ePyramid:
3291  {
3292  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3293  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3294  returnval.push_back(bkey);
3295  returnval.push_back(bkey);
3296 
3297  const LibUtilities::PointsKey pkey1(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha2Beta0);
3298  LibUtilities::BasisKey bkey1(LibUtilities::eModified_C, nummodes, pkey1);
3299  returnval.push_back(bkey1);
3300  }
3301  break;
3302  case LibUtilities::ePrism:
3303  {
3304  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
3305  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
3306  returnval.push_back(bkey);
3307  returnval.push_back(bkey);
3308 
3309  const LibUtilities::PointsKey pkey1(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha1Beta0);
3310  LibUtilities::BasisKey bkey1(LibUtilities::eModified_B, nummodes, pkey1);
3311  returnval.push_back(bkey1);
3312 
3313  }
3314  break;
3315  default:
3316  {
3317  ASSERTL0(false,"Expansion not defined in switch for this shape");
3318  }
3319  break;
3320  }
3321  }
3322  break;
3323 
3324  case eGLL_Lagrange:
3325  {
3326  switch(shape)
3327  {
3328  case LibUtilities::eSegment:
3329  {
3330  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3331  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3332  returnval.push_back(bkey);
3333  }
3334  break;
3335  case LibUtilities::eQuadrilateral:
3336  {
3337  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3338  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3339  returnval.push_back(bkey);
3340  returnval.push_back(bkey);
3341  }
3342  break;
3343  case LibUtilities::eTriangle: // define with corrects points key
3344  // and change to Ortho on construction
3345  {
3346  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3347  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3348  returnval.push_back(bkey);
3349 
3350  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussRadauMAlpha1Beta0);
3351  LibUtilities::BasisKey bkey1(LibUtilities::eOrtho_B, nummodes, pkey1);
3352  returnval.push_back(bkey1);
3353  }
3354  break;
3355  case LibUtilities::eHexahedron:
3356  {
3357  const LibUtilities::PointsKey pkey(nummodes+1,LibUtilities::eGaussLobattoLegendre);
3358  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3359 
3360  returnval.push_back(bkey);
3361  returnval.push_back(bkey);
3362  returnval.push_back(bkey);
3363  }
3364  break;
3365  default:
3366  {
3367  ASSERTL0(false, "Expansion not defined in switch for this shape");
3368  }
3369  break;
3370  }
3371  }
3372  break;
3373 
3374  case eGauss_Lagrange:
3375  {
3376  switch (shape)
3377  {
3378  case LibUtilities::eSegment:
3379  {
3380  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussLegendre);
3381  LibUtilities::BasisKey bkey(LibUtilities::eGauss_Lagrange, nummodes, pkey);
3382 
3383  returnval.push_back(bkey);
3384  }
3385  break;
3386  case LibUtilities::eQuadrilateral:
3387  {
3388  const LibUtilities::PointsKey pkey(nummodes,LibUtilities::eGaussGaussLegendre);
3389  LibUtilities::BasisKey bkey(LibUtilities::eGauss_Lagrange, nummodes, pkey);
3390 
3391  returnval.push_back(bkey);
3392  returnval.push_back(bkey);
3393  }
3394  break;
3395  case LibUtilities::eHexahedron:
3396  {
3397  const LibUtilities::PointsKey pkey(nummodes,LibUtilities::eGaussGaussLegendre);
3398  LibUtilities::BasisKey bkey(LibUtilities::eGauss_Lagrange, nummodes, pkey);
3399 
3400  returnval.push_back(bkey);
3401  returnval.push_back(bkey);
3402  returnval.push_back(bkey);
3403  }
3404  break;
3405  default:
3406  {
3407  ASSERTL0(false, "Expansion not defined in switch for this shape");
3408  }
3409  break;
3410  }
3411  }
3412  break;
3413 
3414  case eOrthogonal:
3415  {
3416  switch (shape)
3417  {
3418  case LibUtilities::eSegment:
3419  {
3420  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3421  LibUtilities::BasisKey bkey(LibUtilities::eOrtho_A, nummodes, pkey);
3422 
3423  returnval.push_back(bkey);
3424  }
3425  break;
3426  case LibUtilities::eTriangle:
3427  {
3428  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3429  LibUtilities::BasisKey bkey(LibUtilities::eOrtho_A, nummodes, pkey);
3430 
3431  returnval.push_back(bkey);
3432 
3433  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussRadauMAlpha1Beta0);
3434  LibUtilities::BasisKey bkey1(LibUtilities::eOrtho_B, nummodes, pkey1);
3435 
3436  returnval.push_back(bkey1);
3437  }
3438  break;
3439  case LibUtilities::eQuadrilateral:
3440  {
3441  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3442  LibUtilities::BasisKey bkey(LibUtilities::eOrtho_A, nummodes, pkey);
3443 
3444  returnval.push_back(bkey);
3445  returnval.push_back(bkey);
3446  }
3447  break;
3448  case LibUtilities::eTetrahedron:
3449  {
3450  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3451  LibUtilities::BasisKey bkey(LibUtilities::eOrtho_A, nummodes, pkey);
3452 
3453  returnval.push_back(bkey);
3454 
3455  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussRadauMAlpha1Beta0);
3456  LibUtilities::BasisKey bkey1(LibUtilities::eOrtho_B, nummodes, pkey1);
3457 
3458  returnval.push_back(bkey1);
3459 
3460  const LibUtilities::PointsKey pkey2(nummodes, LibUtilities::eGaussRadauMAlpha2Beta0);
3461  LibUtilities::BasisKey bkey2(LibUtilities::eOrtho_C, nummodes, pkey2);
3462  }
3463  break;
3464  default:
3465  {
3466  ASSERTL0(false,"Expansion not defined in switch for this shape");
3467  }
3468  break;
3469  }
3470  }
3471  break;
3472 
3473  case eGLL_Lagrange_SEM:
3474  {
3475  switch (shape)
3476  {
3477  case LibUtilities::eSegment:
3478  {
3479  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussLobattoLegendre);
3480  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3481 
3482  returnval.push_back(bkey);
3483  }
3484  break;
3485  case LibUtilities::eQuadrilateral:
3486  {
3487  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussLobattoLegendre);
3488  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3489 
3490  returnval.push_back(bkey);
3491  returnval.push_back(bkey);
3492  }
3493  break;
3494  case LibUtilities::eHexahedron:
3495  {
3496  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussLobattoLegendre);
3497  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
3498 
3499  returnval.push_back(bkey);
3500  returnval.push_back(bkey);
3501  returnval.push_back(bkey);
3502  }
3503  break;
3504  default:
3505  {
3506  ASSERTL0(false,"Expansion not defined in switch for this shape");
3507  }
3508  break;
3509  }
3510  }
3511  break;
3512 
3513 
3514  case eFourier:
3515  {
3516  switch (shape)
3517  {
3518  case LibUtilities::eSegment:
3519  {
3520  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3521  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3522  returnval.push_back(bkey);
3523  }
3524  break;
3525  case LibUtilities::eQuadrilateral:
3526  {
3527  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3528  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3529  returnval.push_back(bkey);
3530  returnval.push_back(bkey);
3531  }
3532  break;
3533  case LibUtilities::eHexahedron:
3534  {
3535  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3536  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3537  returnval.push_back(bkey);
3538  returnval.push_back(bkey);
3539  returnval.push_back(bkey);
3540  }
3541  break;
3542  default:
3543  {
3544  ASSERTL0(false,"Expansion not defined in switch for this shape");
3545  }
3546  break;
3547  }
3548  }
3549  break;
3550 
3551 
3552  case eFourierSingleMode:
3553  {
3554  switch (shape)
3555  {
3556  case LibUtilities::eSegment:
3557  {
3558  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3559  LibUtilities::BasisKey bkey(LibUtilities::eFourierSingleMode, nummodes, pkey);
3560  returnval.push_back(bkey);
3561  }
3562  break;
3563  case LibUtilities::eQuadrilateral:
3564  {
3565  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3566  LibUtilities::BasisKey bkey(LibUtilities::eFourierSingleMode, nummodes, pkey);
3567  returnval.push_back(bkey);
3568  returnval.push_back(bkey);
3569  }
3570  break;
3571  case LibUtilities::eHexahedron:
3572  {
3573  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3574  LibUtilities::BasisKey bkey(LibUtilities::eFourierSingleMode, nummodes, pkey);
3575  returnval.push_back(bkey);
3576  returnval.push_back(bkey);
3577  returnval.push_back(bkey);
3578  }
3579  break;
3580  default:
3581  {
3582  ASSERTL0(false,"Expansion not defined in switch for this shape");
3583  }
3584  break;
3585  }
3586  }
3587  break;
3588 
3589  case eFourierHalfModeRe:
3590  {
3591  switch (shape)
3592  {
3593  case LibUtilities::eSegment:
3594  {
3595  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3596  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeRe, nummodes, pkey);
3597  returnval.push_back(bkey);
3598  }
3599  break;
3600  case LibUtilities::eQuadrilateral:
3601  {
3602  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3603  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeRe, nummodes, pkey);
3604  returnval.push_back(bkey);
3605  returnval.push_back(bkey);
3606  }
3607  break;
3608  case LibUtilities::eHexahedron:
3609  {
3610  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3611  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeRe, nummodes, pkey);
3612  returnval.push_back(bkey);
3613  returnval.push_back(bkey);
3614  returnval.push_back(bkey);
3615  }
3616  break;
3617  default:
3618  {
3619  ASSERTL0(false,"Expansion not defined in switch for this shape");
3620  }
3621  break;
3622  }
3623  }
3624  break;
3625 
3626  case eFourierHalfModeIm:
3627  {
3628  switch (shape)
3629  {
3630  case LibUtilities::eSegment:
3631  {
3632  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3633  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeIm, nummodes, pkey);
3634  returnval.push_back(bkey);
3635  }
3636  break;
3637  case LibUtilities::eQuadrilateral:
3638  {
3639  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3640  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeIm, nummodes, pkey);
3641  returnval.push_back(bkey);
3642  returnval.push_back(bkey);
3643  }
3644  break;
3645  case LibUtilities::eHexahedron:
3646  {
3647  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3648  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeIm, nummodes, pkey);
3649  returnval.push_back(bkey);
3650  returnval.push_back(bkey);
3651  returnval.push_back(bkey);
3652  }
3653  break;
3654  default:
3655  {
3656  ASSERTL0(false,"Expansion not defined in switch for this shape");
3657  }
3658  break;
3659  }
3660  }
3661  break;
3662 
3663  case eChebyshev:
3664  {
3665  switch (shape)
3666  {
3667  case LibUtilities::eSegment:
3668  {
3669  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussChebyshev);
3670  LibUtilities::BasisKey bkey(LibUtilities::eChebyshev, nummodes, pkey);
3671  returnval.push_back(bkey);
3672  }
3673  break;
3674  case LibUtilities::eQuadrilateral:
3675  {
3676  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussChebyshev);
3677  LibUtilities::BasisKey bkey(LibUtilities::eChebyshev, nummodes, pkey);
3678  returnval.push_back(bkey);
3679  returnval.push_back(bkey);
3680  }
3681  break;
3682  case LibUtilities::eHexahedron:
3683  {
3684  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussChebyshev);
3685  LibUtilities::BasisKey bkey(LibUtilities::eChebyshev, nummodes, pkey);
3686  returnval.push_back(bkey);
3687  returnval.push_back(bkey);
3688  returnval.push_back(bkey);
3689  }
3690  break;
3691  default:
3692  {
3693  ASSERTL0(false,"Expansion not defined in switch for this shape");
3694  }
3695  break;
3696  }
3697  }
3698  break;
3699 
3700  case eFourierChebyshev:
3701  {
3702  switch (shape)
3703  {
3704  case LibUtilities::eQuadrilateral:
3705  {
3706  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3707  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3708  returnval.push_back(bkey);
3709 
3710  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussGaussChebyshev);
3711  LibUtilities::BasisKey bkey1(LibUtilities::eChebyshev, nummodes, pkey1);
3712  returnval.push_back(bkey1);
3713  }
3714  break;
3715  default:
3716  {
3717  ASSERTL0(false,"Expansion not defined in switch for this shape");
3718  }
3719  break;
3720  }
3721  }
3722  break;
3723 
3724  case eChebyshevFourier:
3725  {
3726  switch (shape)
3727  {
3728  case LibUtilities::eQuadrilateral:
3729  {
3730  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussGaussChebyshev);
3731  LibUtilities::BasisKey bkey1(LibUtilities::eChebyshev, nummodes, pkey1);
3732  returnval.push_back(bkey1);
3733 
3734  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3735  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3736  returnval.push_back(bkey);
3737  }
3738  break;
3739  default:
3740  {
3741  ASSERTL0(false,"Expansion not defined in switch for this shape");
3742  }
3743  break;
3744  }
3745  }
3746  break;
3747 
3748  case eFourierModified:
3749  {
3750  switch (shape)
3751  {
3752  case LibUtilities::eQuadrilateral:
3753  {
3754  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3755  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3756  returnval.push_back(bkey);
3757 
3758  const LibUtilities::PointsKey pkey1(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3759  LibUtilities::BasisKey bkey1(LibUtilities::eModified_A, nummodes, pkey1);
3760  returnval.push_back(bkey1);
3761  }
3762  break;
3763  default:
3764  {
3765  ASSERTL0(false,"Expansion not defined in switch for this shape");
3766  }
3767  break;
3768  }
3769  }
3770  break;
3771 
3772  default:
3773  {
3774  ASSERTL0(false,"Expansion type not defined");
3775  }
3776  break;
3777 
3778  }
3779 
3780  return returnval;
3781  }
3782 
3783  /**
3784  *
3785  */
3786  LibUtilities::BasisKeyVector
3787  MeshGraph::DefineBasisKeyFromExpansionTypeHomo(
3788  GeometrySharedPtr in,
3789  ExpansionType type_x,
3790  ExpansionType type_y,
3791  ExpansionType type_z,
3792  const int nummodes_x,
3793  const int nummodes_y,
3794  const int nummodes_z)
3795  {
3796  LibUtilities::BasisKeyVector returnval;
3797 
3798  LibUtilities::ShapeType shape = in->GetShapeType();
3799 
3800  switch (shape)
3801  {
3802  case LibUtilities::eSegment:
3803  {
3804  ASSERTL0(false,"Homogeneous expansion not defined for this shape");
3805  }
3806  break;
3807 
3808  case LibUtilities::eQuadrilateral:
3809  {
3810  ASSERTL0(false,"Homogeneous expansion not defined for this shape");
3811  }
3812  break;
3813 
3814  case LibUtilities::eHexahedron:
3815  {
3816  switch(type_x)
3817  {
3818  case eFourier:
3819  {
3820  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eFourierEvenlySpaced);
3821  LibUtilities::BasisKey bkey1(LibUtilities::eFourier,nummodes_x,pkey1);
3822  returnval.push_back(bkey1);
3823  }
3824  break;
3825 
3826  case eFourierSingleMode:
3827  {
3828  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eFourierSingleModeSpaced);
3829  LibUtilities::BasisKey bkey1(LibUtilities::eFourierSingleMode,nummodes_x,pkey1);
3830  returnval.push_back(bkey1);
3831  }
3832  break;
3833 
3834  case eFourierHalfModeRe:
3835  {
3836  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eFourierSingleModeSpaced);
3837  LibUtilities::BasisKey bkey1(LibUtilities::eFourierHalfModeRe,nummodes_x,pkey1);
3838  returnval.push_back(bkey1);
3839  }
3840  break;
3841 
3842  case eFourierHalfModeIm:
3843  {
3844  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eFourierSingleModeSpaced);
3845  LibUtilities::BasisKey bkey1(LibUtilities::eFourierHalfModeIm,nummodes_x,pkey1);
3846  returnval.push_back(bkey1);
3847  }
3848  break;
3849 
3850 
3851  case eChebyshev:
3852  {
3853  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eGaussGaussChebyshev);
3854  LibUtilities::BasisKey bkey1(LibUtilities::eChebyshev,nummodes_x,pkey1);
3855  returnval.push_back(bkey1);
3856  }
3857  break;
3858 
3859 
3860 
3861  default:
3862  {
3863  ASSERTL0(false,"Homogeneous expansion can be of Fourier or Chebyshev type only");
3864  }
3865  break;
3866  }
3867 
3868 
3869  switch(type_y)
3870  {
3871  case eFourier:
3872  {
3873  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eFourierEvenlySpaced);
3874  LibUtilities::BasisKey bkey2(LibUtilities::eFourier,nummodes_y,pkey2);
3875  returnval.push_back(bkey2);
3876  }
3877  break;
3878 
3879 
3880  case eFourierSingleMode:
3881  {
3882  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eFourierSingleModeSpaced);
3883  LibUtilities::BasisKey bkey2(LibUtilities::eFourierSingleMode,nummodes_y,pkey2);
3884  returnval.push_back(bkey2);
3885  }
3886  break;
3887 
3888  case eFourierHalfModeRe:
3889  {
3890  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eFourierSingleModeSpaced);
3891  LibUtilities::BasisKey bkey2(LibUtilities::eFourierHalfModeRe,nummodes_y,pkey2);
3892  returnval.push_back(bkey2);
3893  }
3894  break;
3895 
3896  case eFourierHalfModeIm:
3897  {
3898  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eFourierSingleModeSpaced);
3899  LibUtilities::BasisKey bkey2(LibUtilities::eFourierHalfModeIm,nummodes_y,pkey2);
3900  returnval.push_back(bkey2);
3901  }
3902  break;
3903 
3904  case eChebyshev:
3905  {
3906  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eGaussGaussChebyshev);
3907  LibUtilities::BasisKey bkey2(LibUtilities::eChebyshev,nummodes_y,pkey2);
3908  returnval.push_back(bkey2);
3909  }
3910  break;
3911 
3912  default:
3913  {
3914  ASSERTL0(false,"Homogeneous expansion can be of Fourier or Chebyshev type only");
3915  }
3916  break;
3917  }
3918 
3919  switch(type_z)
3920  {
3921  case eFourier:
3922  {
3923  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eFourierEvenlySpaced);
3924  LibUtilities::BasisKey bkey3(LibUtilities::eFourier,nummodes_z,pkey3);
3925  returnval.push_back(bkey3);
3926  }
3927  break;
3928 
3929  case eFourierSingleMode:
3930  {
3931  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eFourierSingleModeSpaced);
3932  LibUtilities::BasisKey bkey3(LibUtilities::eFourierSingleMode,nummodes_z,pkey3);
3933  returnval.push_back(bkey3);
3934  }
3935  break;
3936 
3937  case eFourierHalfModeRe:
3938  {
3939  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eFourierSingleModeSpaced);
3940  LibUtilities::BasisKey bkey3(LibUtilities::eFourierHalfModeRe,nummodes_z,pkey3);
3941  returnval.push_back(bkey3);
3942  }
3943  break;
3944 
3945  case eFourierHalfModeIm:
3946  {
3947  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eFourierSingleModeSpaced);
3948  LibUtilities::BasisKey bkey3(LibUtilities::eFourierHalfModeIm,nummodes_z,pkey3);
3949  returnval.push_back(bkey3);
3950  }
3951  break;
3952 
3953  case eChebyshev:
3954  {
3955  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eGaussGaussChebyshev);
3956  LibUtilities::BasisKey bkey3(LibUtilities::eChebyshev,nummodes_z,pkey3);
3957  returnval.push_back(bkey3);
3958  }
3959  break;
3960 
3961  default:
3962  {
3963  ASSERTL0(false,"Homogeneous expansion can be of Fourier or Chebyshev type only");
3964  }
3965  break;
3966  }
3967  }
3968  break;
3969 
3970  case LibUtilities::eTriangle:
3971  {
3972  ASSERTL0(false,"Homogeneous expansion not defined for this shape");
3973  }
3974  break;
3975 
3976  case LibUtilities::eTetrahedron:
3977  {
3978  ASSERTL0(false,"Homogeneous expansion not defined for this shape");
3979  }
3980  break;
3981 
3982  default:
3983  ASSERTL0(false,"Expansion not defined in switch for this shape");
3984  break;
3985  }
3986 
3987  return returnval;
3988  }
3989 
3990 
3991  /**
3992  *
3993  */
3994  PointGeomSharedPtr MeshGraph::AddVertex(NekDouble x, NekDouble y, NekDouble z)
3995  {
3996  unsigned int nextId = m_vertSet.rbegin()->first + 1;
3997  PointGeomSharedPtr vert(MemoryManager<PointGeom>::AllocateSharedPtr(m_spaceDimension, nextId, x, y, z));
3998  m_vertSet[nextId] = vert;
3999  return vert;
4000  }
4001 
4002  /**
4003  *
4004  */
4005  SegGeomSharedPtr MeshGraph::AddEdge(PointGeomSharedPtr v0, PointGeomSharedPtr v1,
4006  CurveSharedPtr curveDefinition)
4007  {
4008  PointGeomSharedPtr vertices[] = {v0, v1};
4009  SegGeomSharedPtr edge;
4010  int edgeId = m_segGeoms.rbegin()->first + 1;
4011 
4012  if( curveDefinition )
4013  {
4014  edge = MemoryManager<SegGeom>::AllocateSharedPtr(edgeId, m_spaceDimension, vertices, curveDefinition);
4015  }
4016  else
4017  {
4018  edge = MemoryManager<SegGeom>::AllocateSharedPtr(edgeId, m_spaceDimension, vertices);
4019  }
4020  m_segGeoms[edgeId] = edge;
4021  return edge;
4022  }
4023 
4024 
4025  /**
4026  *
4027  */
4028  TriGeomSharedPtr MeshGraph::AddTriangle(SegGeomSharedPtr edges[], StdRegions::Orientation orient[])
4029  {
4030  int indx = m_triGeoms.rbegin()->first + 1;
4031  TriGeomSharedPtr trigeom(MemoryManager<TriGeom>::AllocateSharedPtr(indx, edges, orient));
4032  trigeom->SetGlobalID(indx);
4033 
4034  m_triGeoms[indx] = trigeom;
4035 
4036  return trigeom;
4037  }
4038 
4039 
4040  /**
4041  *
4042  */
4043  QuadGeomSharedPtr MeshGraph::AddQuadrilateral(SegGeomSharedPtr edges[], StdRegions::Orientation orient[])
4044  {
4045  int indx = m_quadGeoms.rbegin()->first + 1;
4046  QuadGeomSharedPtr quadgeom(MemoryManager<QuadGeom>::AllocateSharedPtr(indx, edges, orient));
4047  quadgeom->SetGlobalID(indx);
4048 
4049  m_quadGeoms[indx] = quadgeom;
4050  return quadgeom;
4051  }
4052 
4053 
4054  /**
4055  *
4056  */
4057  PrismGeomSharedPtr MeshGraph::AddPrism(TriGeomSharedPtr tfaces[PrismGeom::kNtfaces],
4058  QuadGeomSharedPtr qfaces[PrismGeom::kNqfaces])
4059  {
4060  // Setting the orientation is disabled in the reader. Why?
4061  Geometry2DSharedPtr faces[] = { qfaces[0], tfaces[0], qfaces[1], tfaces[1], qfaces[2] };
4062  unsigned int index = m_prismGeoms.rbegin()->first + 1;
4063  PrismGeomSharedPtr prismgeom(MemoryManager<PrismGeom>::AllocateSharedPtr(faces));
4064  prismgeom->SetGlobalID(index);
4065 
4066  m_prismGeoms[index] = prismgeom;
4067  return prismgeom;
4068  }
4069 
4070 
4071  /**
4072  *
4073  */
4074  TetGeomSharedPtr MeshGraph::AddTetrahedron(TriGeomSharedPtr tfaces[TetGeom::kNtfaces])
4075  {
4076  unsigned int index = m_tetGeoms.rbegin()->first + 1;
4077  TetGeomSharedPtr tetgeom(MemoryManager<TetGeom>::AllocateSharedPtr(tfaces));
4078  tetgeom->SetGlobalID(index);
4079 
4080  m_tetGeoms[index] = tetgeom;
4081  return tetgeom;
4082  }
4083 
4084 
4085  /**
4086  *
4087  */
4088  PyrGeomSharedPtr MeshGraph::AddPyramid(TriGeomSharedPtr tfaces[PyrGeom::kNtfaces],
4089  QuadGeomSharedPtr qfaces[PyrGeom::kNqfaces])
4090  {
4091  Geometry2DSharedPtr faces[] = { qfaces[0], tfaces[0], tfaces[1], tfaces[2], tfaces[3] };
4092  unsigned int index = m_pyrGeoms.rbegin()->first + 1;
4093 
4094  PyrGeomSharedPtr pyrgeom(MemoryManager<PyrGeom>::AllocateSharedPtr(faces));
4095  pyrgeom->SetGlobalID(index);
4096 
4097  m_pyrGeoms[index] = pyrgeom;
4098  return pyrgeom;
4099  }
4100 
4101 
4102  /**
4103  *
4104  */
4105  HexGeomSharedPtr MeshGraph::AddHexahedron(QuadGeomSharedPtr qfaces[HexGeom::kNqfaces])
4106  {
4107  unsigned int index = m_hexGeoms.rbegin()->first + 1;
4108  HexGeomSharedPtr hexgeom(MemoryManager<HexGeom>::AllocateSharedPtr(qfaces));
4109  hexgeom->SetGlobalID(index);
4110  m_hexGeoms[index] = hexgeom;
4111  return hexgeom;
4112  }
4113 
4114 
4115  /**
4116  * Generate a single vector of Expansion structs mapping global element
4117  * ID to a corresponding Geometry shared pointer and basis key.
4118  *
4119  * Expansion map ensures elements which appear in multiple composites
4120  * within the domain are only listed once.
4121  */
4122  ExpansionMapShPtr MeshGraph::SetUpExpansionMap(void)
4123  {
4124  ExpansionMapShPtr returnval;
4125  returnval = MemoryManager<ExpansionMap>::AllocateSharedPtr();
4126 
4127  for(int d = 0; d < m_domain.size(); ++d)
4128  {
4129  CompositeMap::const_iterator compIter;
4130 
4131  for (compIter = m_domain[d].begin(); compIter != m_domain[d].end(); ++compIter)
4132  {
4133  GeometryVector::const_iterator x;
4134  for (x = compIter->second->begin(); x != compIter->second->end(); ++x)
4135  {
4136  LibUtilities::BasisKeyVector def;
4137  ExpansionShPtr expansionElementShPtr =
4138  MemoryManager<Expansion>::AllocateSharedPtr(*x, def);
4139  int id = (*x)->GetGlobalID();
4140  (*returnval)[id] = expansionElementShPtr;
4141  }
4142  }
4143  }
4144 
4145  return returnval;
4146  }
4147  }; //end of namespace
4148 }; //end of namespace
StdPrismExp.h
Nektar::SpatialDomains::MeshGraph::SetExpansions
void SetExpansions(std::vector< LibUtilities::FieldDefinitionsSharedPtr > &fielddef)
Sets expansions given field definitions.
Definition: MeshGraph.cpp:2342
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:3094
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:4057
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:3994
Nektar::SpatialDomains::MeshGraph::AddQuadrilateral
QuadGeomSharedPtr AddQuadrilateral(SegGeomSharedPtr edges[], StdRegions::Orientation orient[])
Definition: MeshGraph.cpp:4043
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:3195
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:3127
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:4028
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:4005
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:4088
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:3787
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:3174
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:4105
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:4122
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:3052
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:4074
StdTetExp.h
Nektar::SpatialDomains::ExpansionMapIter
std::map< int, ExpansionShPtr >::iterator ExpansionMapIter
Definition: MeshGraph.h:175
Generated on Fri Dec 2 2016 13:10:49 for Nektar++ by   doxygen 1.8.8