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