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  ASSERTL0(err==TIXML_SUCCESS, "Unable to read mesh dimension.");
136  break;
137  }
138  else
139  {
140  std::string errstr("Unknown attribute: ");
141  errstr += attrName;
142  ASSERTL0(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::TINYXML_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::TINYXML_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::TINYXML_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::TINYXML_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::TINYXML_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  m_domainRange->m_checkShape = false;
1384 
1385  if(m_domainRange == NullDomainRangeShPtr)
1386  {
1387  m_domainRange = MemoryManager<DomainRange>::AllocateSharedPtr();
1388  m_domainRange->m_doXrange = true;
1389  }
1390 
1391  m_domainRange->m_xmin = xmin;
1392  m_domainRange->m_xmax = xmax;
1393 
1394  if(ymin == NekConstants::kNekUnsetDouble)
1395  {
1396  m_domainRange->m_doYrange = false;
1397  }
1398  else
1399  {
1400  m_domainRange->m_doYrange = true;
1401  m_domainRange->m_ymin = ymin;
1402  m_domainRange->m_ymax = ymax;
1403  }
1404 
1405  if(zmin == NekConstants::kNekUnsetDouble)
1406  {
1407  m_domainRange->m_doZrange = false;
1408  }
1409  else
1410  {
1411  m_domainRange->m_doZrange = true;
1412  m_domainRange->m_zmin = zmin;
1413  m_domainRange->m_zmax = zmax;
1414  }
1415  }
1416 
1417  bool MeshGraph::CheckRange(Geometry2D &geom)
1418  {
1419  bool returnval = true;
1420 
1421  if(m_domainRange != NullDomainRangeShPtr)
1422  {
1423  int nverts = geom.GetNumVerts();
1424  int coordim = geom.GetCoordim();
1425 
1426  // exclude elements outside x range if all vertices not in region
1427  if(m_domainRange->m_doXrange)
1428  {
1429  int ncnt_low = 0;
1430  int ncnt_up = 0;
1431  for(int i = 0; i < nverts; ++i)
1432  {
1433  NekDouble xval = (*geom.GetVertex(i))[0];
1434  if(xval < m_domainRange->m_xmin)
1435  {
1436  ncnt_low++;
1437  }
1438 
1439  if(xval > m_domainRange->m_xmax)
1440  {
1441  ncnt_up++;
1442  }
1443  }
1444 
1445  // check for all verts to be less or greater than
1446  // range so that if element spans thin range then
1447  // it is still included
1448  if((ncnt_up == nverts)||(ncnt_low == nverts))
1449  {
1450  returnval = false;
1451  }
1452  }
1453 
1454  // exclude elements outside y range if all vertices not in region
1455  if(m_domainRange->m_doYrange)
1456  {
1457  int ncnt_low = 0;
1458  int ncnt_up = 0;
1459  for(int i = 0; i < nverts; ++i)
1460  {
1461  NekDouble yval = (*geom.GetVertex(i))[1];
1462  if(yval < m_domainRange->m_ymin)
1463  {
1464  ncnt_low++;
1465  }
1466 
1467  if(yval > m_domainRange->m_ymax)
1468  {
1469  ncnt_up++;
1470  }
1471  }
1472 
1473  // check for all verts to be less or greater than
1474  // range so that if element spans thin range then
1475  // it is still included
1476  if((ncnt_up == nverts)||(ncnt_low == nverts))
1477  {
1478  returnval = false;
1479  }
1480  }
1481 
1482  if(coordim > 2)
1483  {
1484  // exclude elements outside z range if all vertices not in region
1485  if(m_domainRange->m_doZrange)
1486  {
1487  int ncnt_low = 0;
1488  int ncnt_up = 0;
1489 
1490  for(int i = 0; i < nverts; ++i)
1491  {
1492  NekDouble zval = (*geom.GetVertex(i))[2];
1493 
1494  if(zval < m_domainRange->m_zmin)
1495  {
1496  ncnt_low++;
1497  }
1498 
1499  if(zval > m_domainRange->m_zmax)
1500  {
1501  ncnt_up++;
1502  }
1503  }
1504 
1505  // check for all verts to be less or greater than
1506  // range so that if element spans thin range then
1507  // it is still included
1508  if((ncnt_up == nverts)||(ncnt_low == nverts))
1509  {
1510  returnval = false;
1511  }
1512  }
1513  }
1514  }
1515  return returnval;
1516  }
1517 
1518 
1519  /* Domain checker for 3D geometries */
1520  bool MeshGraph::CheckRange(Geometry3D &geom)
1521  {
1522  bool returnval = true;
1523 
1524  if(m_domainRange != NullDomainRangeShPtr)
1525  {
1526  int nverts = geom.GetNumVerts();
1527 
1528  if(m_domainRange->m_doXrange)
1529  {
1530  int ncnt_low = 0;
1531  int ncnt_up = 0;
1532 
1533  for(int i = 0; i < nverts; ++i)
1534  {
1535  NekDouble xval = (*geom.GetVertex(i))[0];
1536  if(xval < m_domainRange->m_xmin)
1537  {
1538  ncnt_low++;
1539  }
1540 
1541  if(xval > m_domainRange->m_xmax)
1542  {
1543  ncnt_up++;
1544  }
1545  }
1546 
1547  // check for all verts to be less or greater than
1548  // range so that if element spans thin range then
1549  // it is still included
1550  if((ncnt_up == nverts)||(ncnt_low == nverts))
1551  {
1552  returnval = false;
1553  }
1554  }
1555 
1556  if(m_domainRange->m_doYrange)
1557  {
1558  int ncnt_low = 0;
1559  int ncnt_up = 0;
1560  for(int i = 0; i < nverts; ++i)
1561  {
1562  NekDouble yval = (*geom.GetVertex(i))[1];
1563  if(yval < m_domainRange->m_ymin)
1564  {
1565  ncnt_low++;
1566  }
1567 
1568  if(yval > m_domainRange->m_ymax)
1569  {
1570  ncnt_up++;
1571  }
1572  }
1573 
1574  // check for all verts to be less or greater than
1575  // range so that if element spans thin range then
1576  // it is still included
1577  if((ncnt_up == nverts)||(ncnt_low == nverts))
1578  {
1579  returnval = false;
1580  }
1581  }
1582 
1583  if(m_domainRange->m_doZrange)
1584  {
1585  int ncnt_low = 0;
1586  int ncnt_up = 0;
1587  for(int i = 0; i < nverts; ++i)
1588  {
1589  NekDouble zval = (*geom.GetVertex(i))[2];
1590 
1591  if(zval < m_domainRange->m_zmin)
1592  {
1593  ncnt_low++;
1594  }
1595 
1596  if(zval > m_domainRange->m_zmax)
1597  {
1598  ncnt_up++;
1599  }
1600  }
1601 
1602  // check for all verts to be less or greater than
1603  // range so that if element spans thin range then
1604  // it is still included
1605  if((ncnt_up == nverts)||(ncnt_low == nverts))
1606  {
1607  returnval = false;
1608  }
1609  }
1610 
1611  if(m_domainRange->m_checkShape)
1612  {
1613  if(geom.GetShapeType() != m_domainRange->m_shapeType)
1614  {
1615  returnval = false;
1616  }
1617  }
1618 
1619  }
1620 
1621  return returnval;
1622  }
1623 
1624  /**
1625  *
1626  */
1627  GeometrySharedPtr MeshGraph::GetCompositeItem(int whichComposite, int whichItem)
1628  {
1629  GeometrySharedPtr returnval;
1630  bool error = false;
1631 
1632  if (whichComposite >= 0 && whichComposite < int(m_meshComposites.size()))
1633  {
1634  if (whichItem >= 0 && whichItem < int(m_meshComposites[whichComposite]->size()))
1635  {
1636  returnval = m_meshComposites[whichComposite]->at(whichItem);
1637  }
1638  else
1639  {
1640  error = true;
1641  }
1642  }
1643  else
1644  {
1645  error = true;
1646  }
1647 
1648  if (error)
1649  {
1650  std::ostringstream errStream;
1651  errStream << "Unable to access composite item [" << whichComposite << "][" << whichItem << "].";
1652 
1653  std::string testStr = errStream.str();
1654 
1655  NEKERROR(ErrorUtil::efatal, testStr.c_str());
1656  }
1657 
1658  return returnval;
1659  }
1660 
1661 
1662  /**
1663  *
1664  */
1665  void MeshGraph::GetCompositeList(const std::string &compositeStr, CompositeMap &compositeVector) const
1666  {
1667  // Parse the composites into a list.
1668  typedef vector<unsigned int> SeqVector;
1669  SeqVector seqVector;
1670  bool parseGood = ParseUtils::GenerateSeqVector(compositeStr.c_str(), seqVector);
1671 
1672  ASSERTL0(parseGood && !seqVector.empty(), (std::string("Unable to read composite index range: ") + compositeStr).c_str());
1673 
1674  SeqVector addedVector; // Vector of those composites already added to compositeVector;
1675  for (SeqVector::iterator iter = seqVector.begin(); iter != seqVector.end(); ++iter)
1676  {
1677  // Only add a new one if it does not already exist in vector.
1678  // Can't go back and delete with a vector, so prevent it from
1679  // being added in the first place.
1680  if (std::find(addedVector.begin(), addedVector.end(), *iter) == addedVector.end())
1681  {
1682 
1683  // If the composite listed is not found and we are working
1684  // on a partitioned mesh, silently ignore it.
1685  if (m_meshComposites.find(*iter) == m_meshComposites.end()
1686  && m_meshPartitioned)
1687  {
1688  continue;
1689  }
1690 
1691  addedVector.push_back(*iter);
1692  Composite composite = GetComposite(*iter);
1693  CompositeMap::iterator compIter;
1694  if (composite)
1695  {
1696  compositeVector[*iter] = composite;
1697  }
1698  else
1699  {
1700  char str[64];
1701  ::sprintf(str, "%d", *iter);
1702  NEKERROR(ErrorUtil::ewarning, (std::string("Undefined composite: ") + str).c_str());
1703 
1704  }
1705  }
1706  }
1707  }
1708 
1709 
1710  /**
1711  *
1712  */
1713  const ExpansionMap &MeshGraph::GetExpansions(const std::string variable)
1714  {
1715  ExpansionMapShPtr returnval;
1716 
1717  if(m_expansionMapShPtrMap.count(variable))
1718  {
1719  returnval = m_expansionMapShPtrMap.find(variable)->second;
1720  }
1721  else
1722  {
1723  if(m_expansionMapShPtrMap.count("DefaultVar") == 0)
1724  {
1725  NEKERROR(ErrorUtil::efatal, (std::string("Unable to find expansion vector definition for field: ")+variable).c_str());
1726  }
1727  returnval = m_expansionMapShPtrMap.find("DefaultVar")->second;
1728  m_expansionMapShPtrMap[variable] = returnval;
1729 
1730  NEKERROR(ErrorUtil::ewarning, (std::string("Using Default variable expansion definition for field: ")+variable).c_str());
1731  }
1732 
1733  return *returnval;
1734  }
1735 
1736 
1737  /**
1738  *
1739  */
1740  ExpansionShPtr MeshGraph::GetExpansion(GeometrySharedPtr geom, const std::string variable)
1741  {
1742  ExpansionMapIter iter;
1743  ExpansionShPtr returnval;
1744 
1745  ExpansionMapShPtr expansionMap = m_expansionMapShPtrMap.find(variable)->second;
1746 
1747  for (iter = expansionMap->begin(); iter!=expansionMap->end(); ++iter)
1748  {
1749  if ((iter->second)->m_geomShPtr == geom)
1750  {
1751  returnval = iter->second;
1752  break;
1753  }
1754  }
1755  return returnval;
1756  }
1757 
1758 
1759  /**
1760  *
1761  */
1762  void MeshGraph::SetExpansions(
1763  std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
1764  {
1765  int i, j, k, cnt, id;
1766  GeometrySharedPtr geom;
1767 
1768  ExpansionMapShPtr expansionMap;
1769 
1770  // Loop over fields and determine unique fields string and
1771  // declare whole expansion list
1772  for(i = 0; i < fielddef.size(); ++i)
1773  {
1774  for(j = 0; j < fielddef[i]->m_fields.size(); ++j)
1775  {
1776  std::string field = fielddef[i]->m_fields[j];
1777  if(m_expansionMapShPtrMap.count(field) == 0)
1778  {
1779  expansionMap = MemoryManager<ExpansionMap>::AllocateSharedPtr();
1780  m_expansionMapShPtrMap[field] = expansionMap;
1781 
1782  // check to see if DefaultVar also not set and if so assign it to this expansion
1783  if(m_expansionMapShPtrMap.count("DefaultVar") == 0)
1784  {
1785  m_expansionMapShPtrMap["DefaultVar"] = expansionMap;
1786  }
1787  }
1788 
1789  // loop over all elements in partition and set expansion
1790  expansionMap = m_expansionMapShPtrMap.find(field)->second;
1791  LibUtilities::BasisKeyVector def;
1792 
1793  for(int d = 0; d < m_domain.size(); ++d)
1794  {
1795  CompositeMap::const_iterator compIter;
1796 
1797  for (compIter = m_domain[d].begin();
1798  compIter != m_domain[d].end(); ++compIter)
1799  {
1800  GeometryVector::const_iterator x;
1801  for (x = compIter->second->begin();
1802  x != compIter->second->end(); ++x)
1803  {
1804  ExpansionShPtr expansionElementShPtr =
1805  MemoryManager<Expansion>::
1806  AllocateSharedPtr(*x, def);
1807  int id = (*x)->GetGlobalID();
1808  (*expansionMap)[id] = expansionElementShPtr;
1809  }
1810  }
1811  }
1812  }
1813  }
1814 
1815  // loop over all elements find the geometry shared ptr and
1816  // set up basiskey vector
1817  for(i = 0; i < fielddef.size(); ++i)
1818  {
1819  cnt = 0;
1820  std::vector<std::string> fields = fielddef[i]->m_fields;
1821  std::vector<unsigned int> nmodes = fielddef[i]->m_numModes;
1822  std::vector<LibUtilities::BasisType> basis = fielddef[i]->m_basis;
1823  bool pointDef = fielddef[i]->m_pointsDef;
1824  bool numPointDef = fielddef[i]->m_numPointsDef;
1825 
1826  // Check points and numpoints
1827  std::vector<unsigned int> npoints = fielddef[i]->m_numPoints;
1828  std::vector<LibUtilities::PointsType> points = fielddef[i]->m_points;
1829 
1830  bool UniOrder = fielddef[i]->m_uniOrder;
1831 
1832  int check = 0;
1833  for (j=0; j< basis.size(); ++j)
1834  {
1835  if ( (strcmp(LibUtilities::BasisTypeMap[basis[j]], "Modified_A") == 0) ||
1836  (strcmp(LibUtilities::BasisTypeMap[basis[j]], "Modified_B") == 0) ||
1837  (strcmp(LibUtilities::BasisTypeMap[basis[j]], "Modified_C") == 0) ||
1838  (strcmp(LibUtilities::BasisTypeMap[basis[j]], "GLL_Lagrange") == 0) ||
1839  (strcmp(LibUtilities::BasisTypeMap[basis[j]], "Gauss_Lagrange") == 0) ||
1840  (strcmp(LibUtilities::BasisTypeMap[basis[j]], "Fourier") == 0) ||
1841  (strcmp(LibUtilities::BasisTypeMap[basis[j]], "FourierSingleMode") == 0)||
1842  (strcmp(LibUtilities::BasisTypeMap[basis[j]], "FourierHalfModeRe") == 0) ||
1843  (strcmp(LibUtilities::BasisTypeMap[basis[j]], "FourierHalfModeIm") == 0))
1844  {
1845  check++;
1846  }
1847  }
1848 
1849  if (check==basis.size())
1850  {
1851  for (j = 0; j < fielddef[i]->m_elementIDs.size(); ++j)
1852  {
1853 
1854  LibUtilities::BasisKeyVector bkeyvec;
1855  id = fielddef[i]->m_elementIDs[j];
1856 
1857  switch (fielddef[i]->m_shapeType)
1858  {
1859  case LibUtilities::eSegment:
1860  {
1861  if(m_segGeoms.count(fielddef[i]->m_elementIDs[j]) == 0)
1862  {
1863  // skip element likely from parallel read
1864  continue;
1865  }
1866  geom = m_segGeoms[fielddef[i]->m_elementIDs[j]];
1867 
1868  LibUtilities::PointsKey pkey(nmodes[cnt]+1, LibUtilities::eGaussLobattoLegendre);
1869 
1870  if(numPointDef&&pointDef)
1871  {
1872  const LibUtilities::PointsKey pkey1(npoints[cnt], points[0]);
1873  pkey = pkey1;
1874  }
1875  else if(!numPointDef&&pointDef)
1876  {
1877  const LibUtilities::PointsKey pkey1(nmodes[cnt]+1, points[0]);
1878  pkey = pkey1;
1879  }
1880  else if(numPointDef&&!pointDef)
1881  {
1882  const LibUtilities::PointsKey pkey1(npoints[cnt], LibUtilities::eGaussLobattoLegendre);
1883  pkey = pkey1;
1884  }
1885 
1886  LibUtilities::BasisKey bkey(basis[0], nmodes[cnt], pkey);
1887 
1888  if(!UniOrder)
1889  {
1890  cnt++;
1891  }
1892  bkeyvec.push_back(bkey);
1893  }
1894  break;
1895  case LibUtilities::eTriangle:
1896  {
1897  if(m_triGeoms.count(fielddef[i]->m_elementIDs[j]) == 0)
1898  {
1899  // skip element likely from parallel read
1900  continue;
1901  }
1902  geom = m_triGeoms[fielddef[i]->m_elementIDs[j]];
1903 
1904  LibUtilities::PointsKey pkey(nmodes[cnt]+1, LibUtilities::eGaussLobattoLegendre);
1905  if(numPointDef&&pointDef)
1906  {
1907  const LibUtilities::PointsKey pkey2(npoints[cnt], points[0]);
1908  pkey = pkey2;
1909  }
1910  else if(!numPointDef&&pointDef)
1911  {
1912  const LibUtilities::PointsKey pkey2(nmodes[cnt]+1, points[0]);
1913  pkey = pkey2;
1914  }
1915  else if(numPointDef&&!pointDef)
1916  {
1917  const LibUtilities::PointsKey pkey2(npoints[cnt], LibUtilities::eGaussLobattoLegendre);
1918  pkey = pkey2;
1919  }
1920  LibUtilities::BasisKey bkey(basis[0], nmodes[cnt], pkey);
1921 
1922  bkeyvec.push_back(bkey);
1923 
1924  LibUtilities::PointsKey pkey1(nmodes[cnt+1], LibUtilities::eGaussRadauMAlpha1Beta0);
1925  if(numPointDef&&pointDef)
1926  {
1927  const LibUtilities::PointsKey pkey2(npoints[cnt+1], points[1]);
1928  pkey1 = pkey2;
1929  }
1930  else if(!numPointDef&&pointDef)
1931  {
1932  const LibUtilities::PointsKey pkey2(nmodes[cnt+1], points[1]);
1933  pkey1 = pkey2;
1934  }
1935  else if(numPointDef&&!pointDef)
1936  {
1937  const LibUtilities::PointsKey pkey2(npoints[cnt+1], LibUtilities::eGaussRadauMAlpha1Beta0);
1938  pkey1 = pkey2;
1939  }
1940  LibUtilities::BasisKey bkey1(basis[1], nmodes[cnt+1], pkey1);
1941  bkeyvec.push_back(bkey1);
1942 
1943  if(!UniOrder)
1944  {
1945  cnt += 2;
1946  }
1947  }
1948  break;
1949  case LibUtilities::eQuadrilateral:
1950  {
1951  if(m_quadGeoms.count(fielddef[i]->m_elementIDs[j]) == 0)
1952  {
1953  // skip element likely from parallel read
1954  continue;
1955  }
1956 
1957  geom = m_quadGeoms[fielddef[i]->m_elementIDs[j]];
1958 
1959  for(int b = 0; b < 2; ++b)
1960  {
1961  LibUtilities::PointsKey pkey(nmodes[cnt+b]+1, LibUtilities::eGaussLobattoLegendre);
1962 
1963  if(numPointDef&&pointDef)
1964  {
1965  const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
1966  pkey = pkey2;
1967  }
1968  else if(!numPointDef&&pointDef)
1969  {
1970  const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
1971  pkey = pkey2;
1972  }
1973  else if(numPointDef&&!pointDef)
1974  {
1975  const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
1976  pkey = pkey2;
1977  }
1978  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
1979  bkeyvec.push_back(bkey);
1980  }
1981 
1982  if(!UniOrder)
1983  {
1984  cnt += 2;
1985  }
1986  }
1987  break;
1988 
1989  case LibUtilities::eTetrahedron:
1990  {
1991  k = fielddef[i]->m_elementIDs[j];
1992 
1993  // allow for possibility that fielddef is
1994  // larger than m_graph which can happen in
1995  // parallel runs
1996  if(m_tetGeoms.count(k) == 0)
1997  {
1998  continue;
1999  }
2000  geom = m_tetGeoms[k];
2001 
2002 #if 0 //all gll
2003  for(int b = 0; b < 3; ++b)
2004  {
2005  LibUtilities::PointsKey pkey(nmodes[cnt+b], LibUtilities::eGaussLobattoLegendre);
2006 
2007  if(numPointDef&&pointDef)
2008  {
2009  const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
2010  pkey = pkey2;
2011  }
2012  else if(!numPointDef&&pointDef)
2013  {
2014  const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
2015  pkey = pkey2;
2016  }
2017  else if(numPointDef&&!pointDef)
2018  {
2019  const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
2020  pkey = pkey2;
2021  }
2022 
2023  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2024 
2025  bkeyvec.push_back(bkey);
2026  }
2027 #else
2028  {
2029  LibUtilities::PointsKey pkey(nmodes[cnt], LibUtilities::eGaussLobattoLegendre);
2030 
2031  if(numPointDef&&pointDef)
2032  {
2033  const LibUtilities::PointsKey pkey2(npoints[cnt],points[0]);
2034  pkey = pkey2;
2035  }
2036  else if(!numPointDef&&pointDef)
2037  {
2038  const LibUtilities::PointsKey pkey2(nmodes[cnt]+1,points[0]);
2039  pkey = pkey2;
2040  }
2041  else if(numPointDef&&!pointDef)
2042  {
2043  const LibUtilities::PointsKey pkey2(npoints[cnt],LibUtilities::eGaussLobattoLegendre);
2044  pkey = pkey2;
2045  }
2046 
2047  LibUtilities::BasisKey bkey(basis[0],nmodes[cnt],pkey);
2048 
2049  bkeyvec.push_back(bkey);
2050  }
2051  {
2052  LibUtilities::PointsKey pkey(nmodes[cnt+1], LibUtilities::eGaussRadauMAlpha1Beta0);
2053 
2054  if(numPointDef&&pointDef)
2055  {
2056  const LibUtilities::PointsKey pkey2(npoints[cnt+1],points[1]);
2057  pkey = pkey2;
2058  }
2059  else if(!numPointDef&&pointDef)
2060  {
2061  const LibUtilities::PointsKey pkey2(nmodes[cnt+1]+1,points[1]);
2062  pkey = pkey2;
2063  }
2064  else if(numPointDef&&!pointDef)
2065  {
2066  const LibUtilities::PointsKey pkey2(npoints[cnt+1],LibUtilities::eGaussRadauMAlpha1Beta0);
2067  pkey = pkey2;
2068  }
2069 
2070  LibUtilities::BasisKey bkey(basis[1],nmodes[cnt+1],pkey);
2071 
2072  bkeyvec.push_back(bkey);
2073  }
2074 
2075  {
2076  LibUtilities::PointsKey pkey(nmodes[cnt+2], LibUtilities::eGaussRadauMAlpha2Beta0);
2077 
2078  if(numPointDef&&pointDef)
2079  {
2080  const LibUtilities::PointsKey pkey2(npoints[cnt+2],points[2]);
2081  pkey = pkey2;
2082  }
2083  else if(!numPointDef&&pointDef)
2084  {
2085  const LibUtilities::PointsKey pkey2(nmodes[cnt+2]+1,points[2]);
2086  pkey = pkey2;
2087  }
2088  else if(numPointDef&&!pointDef)
2089  {
2090  const LibUtilities::PointsKey pkey2(npoints[cnt+2],LibUtilities::eGaussRadauMAlpha1Beta0);
2091  pkey = pkey2;
2092  }
2093 
2094  LibUtilities::BasisKey bkey(basis[2],nmodes[cnt+2],pkey);
2095 
2096  bkeyvec.push_back(bkey);
2097  }
2098 #endif
2099 
2100  if(!UniOrder)
2101  {
2102  cnt += 3;
2103  }
2104  }
2105  break;
2106  case LibUtilities::ePrism:
2107  {
2108  k = fielddef[i]->m_elementIDs[j];
2109  if(m_prismGeoms.count(k) == 0)
2110  {
2111  continue;
2112  }
2113  geom = m_prismGeoms[k];
2114 
2115 #if 0 // all GLL
2116  for(int b = 0; b < 3; ++b)
2117  {
2118  LibUtilities::PointsKey pkey(nmodes[cnt+b],LibUtilities::eGaussLobattoLegendre);
2119 
2120  if(numPointDef&&pointDef)
2121  {
2122  const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
2123  pkey = pkey2;
2124  }
2125  else if(!numPointDef&&pointDef)
2126  {
2127  const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
2128  pkey = pkey2;
2129  }
2130  else if(numPointDef&&!pointDef)
2131  {
2132  const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
2133  pkey = pkey2;
2134  }
2135 
2136  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2137  bkeyvec.push_back(bkey);
2138  }
2139 #else
2140  for(int b = 0; b < 2; ++b)
2141  {
2142  LibUtilities::PointsKey pkey(nmodes[cnt+b],LibUtilities::eGaussLobattoLegendre);
2143 
2144  if(numPointDef&&pointDef)
2145  {
2146  const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
2147  pkey = pkey2;
2148  }
2149  else if(!numPointDef&&pointDef)
2150  {
2151  const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
2152  pkey = pkey2;
2153  }
2154  else if(numPointDef&&!pointDef)
2155  {
2156  const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
2157  pkey = pkey2;
2158  }
2159 
2160  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2161  bkeyvec.push_back(bkey);
2162  }
2163 
2164  {
2165  LibUtilities::PointsKey pkey(nmodes[cnt+2],LibUtilities::eGaussRadauMAlpha1Beta0);
2166 
2167  if(numPointDef&&pointDef)
2168  {
2169  const LibUtilities::PointsKey pkey2(npoints[cnt+2],points[2]);
2170  pkey = pkey2;
2171  }
2172  else if(!numPointDef&&pointDef)
2173  {
2174  const LibUtilities::PointsKey pkey2(nmodes[cnt+2]+1,points[2]);
2175  pkey = pkey2;
2176  }
2177  else if(numPointDef&&!pointDef)
2178  {
2179  const LibUtilities::PointsKey pkey2(npoints[cnt+2],LibUtilities::eGaussLobattoLegendre);
2180  pkey = pkey2;
2181  }
2182 
2183  LibUtilities::BasisKey bkey(basis[2],nmodes[cnt+2],pkey);
2184  bkeyvec.push_back(bkey);
2185  }
2186 #endif
2187  if(!UniOrder)
2188  {
2189  cnt += 3;
2190  }
2191  }
2192  break;
2193  case LibUtilities::ePyramid:
2194  {
2195  k = fielddef[i]->m_elementIDs[j];
2196  ASSERTL0(m_pyrGeoms.find(k) != m_pyrGeoms.end(),
2197  "Failed to find geometry with same global id");
2198  geom = m_pyrGeoms[k];
2199 
2200  for(int b = 0; b < 3; ++b)
2201  {
2202  LibUtilities::PointsKey pkey(nmodes[cnt+b],points[b]);
2203 
2204  if(numPointDef&&pointDef)
2205  {
2206  const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
2207  pkey = pkey2;
2208  }
2209  else if(!numPointDef&&pointDef)
2210  {
2211  const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
2212  pkey = pkey2;
2213  }
2214  else if(numPointDef&&!pointDef)
2215  {
2216  const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
2217  pkey = pkey2;
2218  }
2219 
2220  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2221  bkeyvec.push_back(bkey);
2222  }
2223 
2224  if(!UniOrder)
2225  {
2226  cnt += 3;
2227  }
2228  }
2229  break;
2230  case LibUtilities::eHexahedron:
2231  {
2232  k = fielddef[i]->m_elementIDs[j];
2233  if(m_hexGeoms.count(k) == 0)
2234  {
2235  continue;
2236  }
2237 
2238  geom = m_hexGeoms[k];
2239 
2240  for(int b = 0; b < 3; ++b)
2241  {
2242  LibUtilities::PointsKey pkey(nmodes[cnt+b],LibUtilities::eGaussLobattoLegendre);
2243 
2244  if(numPointDef&&pointDef)
2245  {
2246  const LibUtilities::PointsKey pkey2(npoints[cnt+b],points[b]);
2247  pkey = pkey2;
2248  }
2249  else if(!numPointDef&&pointDef)
2250  {
2251  const LibUtilities::PointsKey pkey2(nmodes[cnt+b]+1,points[b]);
2252  pkey = pkey2;
2253  }
2254  else if(numPointDef&&!pointDef)
2255  {
2256  const LibUtilities::PointsKey pkey2(npoints[cnt+b],LibUtilities::eGaussLobattoLegendre);
2257  pkey = pkey2;
2258  }
2259 
2260  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2261  bkeyvec.push_back(bkey);
2262  }
2263 
2264  if(!UniOrder)
2265  {
2266  cnt += 3;
2267  }
2268  }
2269  break;
2270  default:
2271  ASSERTL0(false,"Need to set up for pyramid and prism 3D Expansions");
2272  break;
2273  }
2274 
2275  for(k = 0; k < fields.size(); ++k)
2276  {
2277  expansionMap = m_expansionMapShPtrMap.find(fields[k])->second;
2278  if((*expansionMap).find(id) != (*expansionMap).end())
2279  {
2280  (*expansionMap)[id]->m_geomShPtr = geom;
2281  (*expansionMap)[id]->m_basisKeyVector = bkeyvec;
2282  }
2283  }
2284  }
2285  }
2286  else
2287  {
2288  ASSERTL0(false,"Need to set up for non Modified basis");
2289  }
2290  }
2291  }
2292 
2293 
2294  /**
2295  *
2296  */
2297  void MeshGraph::SetExpansions(
2298  std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef,
2299  std::vector< std::vector<LibUtilities::PointsType> > &pointstype)
2300  {
2301  int i,j,k,g,h,cnt,id;
2302  GeometrySharedPtr geom;
2303 
2304  ExpansionMapShPtr expansionMap;
2305 
2306  // Loop over fields and determine unique fields string and
2307  // declare whole expansion list
2308  for(i = 0; i < fielddef.size(); ++i)
2309  {
2310  for(j = 0; j < fielddef[i]->m_fields.size(); ++j)
2311  {
2312  std::string field = fielddef[i]->m_fields[j];
2313  if(m_expansionMapShPtrMap.count(field) == 0)
2314  {
2315  expansionMap = MemoryManager<ExpansionMap>::AllocateSharedPtr();
2316  m_expansionMapShPtrMap[field] = expansionMap;
2317 
2318  // check to see if DefaultVar also not set and if so assign it to this expansion
2319  if(m_expansionMapShPtrMap.count("DefaultVar") == 0)
2320  {
2321  m_expansionMapShPtrMap["DefaultVar"] = expansionMap;
2322  }
2323 
2324  // loop over all elements and set expansion
2325  for(k = 0; k < fielddef.size(); ++k)
2326  {
2327  for(h = 0; h < fielddef[k]->m_fields.size(); ++h)
2328  {
2329  if(fielddef[k]->m_fields[h] == field)
2330  {
2331  expansionMap = m_expansionMapShPtrMap.find(field)->second;
2332  LibUtilities::BasisKeyVector def;
2333 
2334  for(g = 0; g < fielddef[k]->m_elementIDs.size(); ++g)
2335  {
2336  ExpansionShPtr tmpexp =
2337  MemoryManager<Expansion>::AllocateSharedPtr(geom, def);
2338  (*expansionMap)[fielddef[k]->m_elementIDs[g]] = tmpexp;
2339  }
2340  }
2341  }
2342  }
2343  }
2344  }
2345  }
2346 
2347 
2348  // loop over all elements find the geometry shared ptr and
2349  // set up basiskey vector
2350  for(i = 0; i < fielddef.size(); ++i)
2351  {
2352  cnt = 0;
2353  std::vector<std::string> fields = fielddef[i]->m_fields;
2354  std::vector<unsigned int> nmodes = fielddef[i]->m_numModes;
2355  std::vector<LibUtilities::BasisType> basis = fielddef[i]->m_basis;
2356  bool UniOrder = fielddef[i]->m_uniOrder;
2357 
2358  for(j = 0; j < fielddef[i]->m_elementIDs.size(); ++j)
2359  {
2360 
2361  LibUtilities::BasisKeyVector bkeyvec;
2362  id = fielddef[i]->m_elementIDs[j];
2363 
2364  switch(fielddef[i]->m_shapeType)
2365  {
2366  case LibUtilities::eSegment:
2367  {
2368  k = fielddef[i]->m_elementIDs[j];
2369  ASSERTL0(m_segGeoms.find(k) != m_segGeoms.end(),
2370  "Failed to find geometry with same global id.");
2371  geom = m_segGeoms[k];
2372 
2373  const LibUtilities::PointsKey pkey(nmodes[cnt], pointstype[i][0]);
2374  LibUtilities::BasisKey bkey(basis[0], nmodes[cnt], pkey);
2375  if(!UniOrder)
2376  {
2377  cnt++;
2378  }
2379  bkeyvec.push_back(bkey);
2380  }
2381  break;
2382  case LibUtilities::eTriangle:
2383  {
2384  k = fielddef[i]->m_elementIDs[j];
2385  ASSERTL0(m_triGeoms.find(k) != m_triGeoms.end(),
2386  "Failed to find geometry with same global id.");
2387  geom = m_triGeoms[k];
2388  for(int b = 0; b < 2; ++b)
2389  {
2390  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
2391  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2392  bkeyvec.push_back(bkey);
2393  }
2394 
2395  if(!UniOrder)
2396  {
2397  cnt += 2;
2398  }
2399  }
2400  break;
2401  case LibUtilities::eQuadrilateral:
2402  {
2403  k = fielddef[i]->m_elementIDs[j];
2404  ASSERTL0(m_quadGeoms.find(k) != m_quadGeoms.end(),
2405  "Failed to find geometry with same global id");
2406  geom = m_quadGeoms[k];
2407 
2408  for(int b = 0; b < 2; ++b)
2409  {
2410  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
2411  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2412  bkeyvec.push_back(bkey);
2413  }
2414 
2415  if(!UniOrder)
2416  {
2417  cnt += 2;
2418  }
2419  }
2420  break;
2421  case LibUtilities::eTetrahedron:
2422  {
2423  k = fielddef[i]->m_elementIDs[j];
2424  ASSERTL0(m_tetGeoms.find(k) != m_tetGeoms.end(),
2425  "Failed to find geometry with same global id");
2426  geom = m_tetGeoms[k];
2427 
2428  for(int b = 0; b < 3; ++b)
2429  {
2430  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
2431  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2432  bkeyvec.push_back(bkey);
2433  }
2434 
2435  if(!UniOrder)
2436  {
2437  cnt += 2;
2438  }
2439  }
2440  break;
2441  case LibUtilities::ePyramid:
2442  {
2443  k = fielddef[i]->m_elementIDs[j];
2444  ASSERTL0(m_pyrGeoms.find(k) != m_pyrGeoms.end(),
2445  "Failed to find geometry with same global id");
2446  geom = m_pyrGeoms[k];
2447 
2448  for(int b = 0; b < 3; ++b)
2449  {
2450  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
2451  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2452  bkeyvec.push_back(bkey);
2453  }
2454 
2455  if(!UniOrder)
2456  {
2457  cnt += 2;
2458  }
2459  }
2460  break;
2461  case LibUtilities::ePrism:
2462  {
2463  k = fielddef[i]->m_elementIDs[j];
2464  ASSERTL0(m_prismGeoms.find(k) != m_prismGeoms.end(),
2465  "Failed to find geometry with same global id");
2466  geom = m_prismGeoms[k];
2467 
2468  for(int b = 0; b < 3; ++b)
2469  {
2470  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
2471  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2472  bkeyvec.push_back(bkey);
2473  }
2474 
2475  if(!UniOrder)
2476  {
2477  cnt += 2;
2478  }
2479  }
2480  break;
2481  case LibUtilities::eHexahedron:
2482  {
2483  k = fielddef[i]->m_elementIDs[j];
2484  ASSERTL0(m_hexGeoms.find(k) != m_hexGeoms.end(),
2485  "Failed to find geometry with same global id");
2486  geom = m_hexGeoms[k];
2487 
2488  for(int b = 0; b < 3; ++b)
2489  {
2490  const LibUtilities::PointsKey pkey(nmodes[cnt+b],pointstype[i][b]);
2491  LibUtilities::BasisKey bkey(basis[b],nmodes[cnt+b],pkey);
2492  bkeyvec.push_back(bkey);
2493  }
2494 
2495  if(!UniOrder)
2496  {
2497  cnt += 2;
2498  }
2499  }
2500  break;
2501  default:
2502  ASSERTL0(false,"Need to set up for pyramid and prism 3D Expansions");
2503  break;
2504  }
2505 
2506  for(k = 0; k < fields.size(); ++k)
2507  {
2508  expansionMap = m_expansionMapShPtrMap.find(fields[k])->second;
2509  if((*expansionMap).find(id) != (*expansionMap).end())
2510  {
2511  (*expansionMap)[id]->m_geomShPtr = geom;
2512  (*expansionMap)[id]->m_basisKeyVector = bkeyvec;
2513  }
2514  }
2515  }
2516  }
2517  }
2518 
2519  /**
2520  * \brief Reset all points keys to have equispaced points with
2521  * optional arguemtn of \a npoints which redefines how many
2522  * points are to be used.
2523  */
2524  void MeshGraph::SetExpansionsToEvenlySpacedPoints(int npoints)
2525  {
2526  ExpansionMapShPtrMapIter it;
2527 
2528  // iterate over all defined expansions
2529  for(it = m_expansionMapShPtrMap.begin(); it != m_expansionMapShPtrMap.end(); ++it)
2530  {
2531  ExpansionMapIter expIt;
2532 
2533  for(expIt = it->second->begin(); expIt != it->second->end(); ++expIt)
2534  {
2535  for(int i = 0; i < expIt->second->m_basisKeyVector.size(); ++i)
2536  {
2537  LibUtilities::BasisKey bkeyold = expIt->second->m_basisKeyVector[i];
2538 
2539  int npts;
2540 
2541  if(npoints) // use input
2542  {
2543  npts = npoints;
2544  }
2545  else
2546  {
2547  npts = bkeyold.GetNumModes();
2548  }
2549 
2550 
2551  const LibUtilities::PointsKey pkey(npts,LibUtilities::ePolyEvenlySpaced);
2552  LibUtilities::BasisKey bkeynew(bkeyold.GetBasisType(),bkeyold.GetNumModes(), pkey);
2553  expIt->second->m_basisKeyVector[i] = bkeynew;
2554 
2555  }
2556  }
2557  }
2558  }
2559 
2560 
2561 
2562  /**
2563  * For each element of shape given by \a shape in field \a
2564  * var, replace the current BasisKeyVector describing the
2565  * expansion in each dimension, with the one provided by \a
2566  * keys.
2567  *
2568  * @TODO: Allow selection of elements through a CompositeVector,
2569  * as well as by type.
2570  *
2571  * @param shape The shape of elements to be changed.
2572  * @param keys The new basis vector to apply to those elements.
2573  */
2574  void MeshGraph::SetBasisKey(LibUtilities::ShapeType shape,
2575  LibUtilities::BasisKeyVector &keys,
2576  std::string var)
2577  {
2578  ExpansionMapIter elemIter;
2579 
2580  ExpansionMapShPtr expansionMap = m_expansionMapShPtrMap.find(var)->second;
2581 
2582  for (elemIter = expansionMap->begin(); elemIter != expansionMap->end(); ++elemIter)
2583  {
2584  if ((elemIter->second)->m_geomShPtr->GetShapeType() == shape)
2585  {
2586  (elemIter->second)->m_basisKeyVector = keys;
2587  }
2588  }
2589  }
2590 
2591 
2592  /**
2593  *
2594  */
2595  LibUtilities::BasisKeyVector MeshGraph::DefineBasisKeyFromExpansionType(
2596  GeometrySharedPtr in,
2597  ExpansionType type,
2598  const int nummodes)
2599  {
2600  LibUtilities::BasisKeyVector returnval;
2601 
2602  LibUtilities::ShapeType shape= in->GetShapeType();
2603 
2604  int quadoffset = 1;
2605  switch(type)
2606  {
2607  case eModified:
2608  quadoffset = 1;
2609  break;
2610  case eModifiedQuadPlus1:
2611  quadoffset = 2;
2612  break;
2613  case eModifiedQuadPlus2:
2614  quadoffset = 3;
2615  break;
2616  default:
2617  break;
2618  }
2619 
2620  switch(type)
2621  {
2622  case eModified:
2623  case eModifiedQuadPlus1:
2624  case eModifiedQuadPlus2:
2625  {
2626  switch (shape)
2627  {
2628  case LibUtilities::eSegment:
2629  {
2630  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
2631  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
2632  returnval.push_back(bkey);
2633  }
2634  break;
2635  case LibUtilities::eQuadrilateral:
2636  {
2637  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
2638  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
2639  returnval.push_back(bkey);
2640  returnval.push_back(bkey);
2641  }
2642  break;
2643  case LibUtilities::eHexahedron:
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  returnval.push_back(bkey);
2649  returnval.push_back(bkey);
2650  }
2651  break;
2652  case LibUtilities::eTriangle:
2653  {
2654  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
2655  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
2656  returnval.push_back(bkey);
2657 
2658  const LibUtilities::PointsKey pkey1(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha1Beta0);
2659  LibUtilities::BasisKey bkey1(LibUtilities::eModified_B, nummodes, pkey1);
2660 
2661  returnval.push_back(bkey1);
2662  }
2663  break;
2664  case LibUtilities::eTetrahedron:
2665  {
2666  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
2667  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
2668  returnval.push_back(bkey);
2669 
2670  const LibUtilities::PointsKey pkey1(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha1Beta0);
2671  LibUtilities::BasisKey bkey1(LibUtilities::eModified_B, nummodes, pkey1);
2672  returnval.push_back(bkey1);
2673 
2674  const LibUtilities::PointsKey pkey2(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha2Beta0);
2675  LibUtilities::BasisKey bkey2(LibUtilities::eModified_C, nummodes, pkey2);
2676  returnval.push_back(bkey2);
2677  }
2678  break;
2679  case LibUtilities::ePyramid:
2680  {
2681  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
2682  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
2683  returnval.push_back(bkey);
2684  returnval.push_back(bkey);
2685 
2686  const LibUtilities::PointsKey pkey1(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha2Beta0);
2687  LibUtilities::BasisKey bkey1(LibUtilities::eModified_C, nummodes, pkey1);
2688  returnval.push_back(bkey1);
2689  }
2690  break;
2691  case LibUtilities::ePrism:
2692  {
2693  const LibUtilities::PointsKey pkey(nummodes+quadoffset, LibUtilities::eGaussLobattoLegendre);
2694  LibUtilities::BasisKey bkey(LibUtilities::eModified_A, nummodes, pkey);
2695  returnval.push_back(bkey);
2696  returnval.push_back(bkey);
2697 
2698  const LibUtilities::PointsKey pkey1(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha1Beta0);
2699  LibUtilities::BasisKey bkey1(LibUtilities::eModified_B, nummodes, pkey1);
2700  returnval.push_back(bkey1);
2701 
2702  }
2703  break;
2704  default:
2705  {
2706  ASSERTL0(false,"Expansion not defined in switch for this shape");
2707  }
2708  break;
2709  }
2710  }
2711  break;
2712 
2713  case eGLL_Lagrange:
2714  {
2715  switch(shape)
2716  {
2717  case LibUtilities::eSegment:
2718  {
2719  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
2720  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
2721  returnval.push_back(bkey);
2722  }
2723  break;
2724  case LibUtilities::eQuadrilateral:
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  returnval.push_back(bkey);
2730  }
2731  break;
2732  case LibUtilities::eTriangle: // define with corrects points key
2733  // and change to Ortho on construction
2734  {
2735  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
2736  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
2737  returnval.push_back(bkey);
2738 
2739  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussRadauMAlpha1Beta0);
2740  LibUtilities::BasisKey bkey1(LibUtilities::eOrtho_B, nummodes, pkey1);
2741  returnval.push_back(bkey1);
2742  }
2743  break;
2744  case LibUtilities::eHexahedron:
2745  {
2746  const LibUtilities::PointsKey pkey(nummodes+1,LibUtilities::eGaussLobattoLegendre);
2747  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
2748 
2749  returnval.push_back(bkey);
2750  returnval.push_back(bkey);
2751  returnval.push_back(bkey);
2752  }
2753  break;
2754  default:
2755  {
2756  ASSERTL0(false, "Expansion not defined in switch for this shape");
2757  }
2758  break;
2759  }
2760  }
2761  break;
2762 
2763  case eGauss_Lagrange:
2764  {
2765  switch (shape)
2766  {
2767  case LibUtilities::eSegment:
2768  {
2769  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussLegendre);
2770  LibUtilities::BasisKey bkey(LibUtilities::eGauss_Lagrange, nummodes, pkey);
2771 
2772  returnval.push_back(bkey);
2773  }
2774  break;
2775  case LibUtilities::eQuadrilateral:
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  }
2783  break;
2784  case LibUtilities::eHexahedron:
2785  {
2786  const LibUtilities::PointsKey pkey(nummodes,LibUtilities::eGaussGaussLegendre);
2787  LibUtilities::BasisKey bkey(LibUtilities::eGauss_Lagrange, nummodes, pkey);
2788 
2789  returnval.push_back(bkey);
2790  returnval.push_back(bkey);
2791  returnval.push_back(bkey);
2792  }
2793  break;
2794  default:
2795  {
2796  ASSERTL0(false, "Expansion not defined in switch for this shape");
2797  }
2798  break;
2799  }
2800  }
2801  break;
2802 
2803  case eOrthogonal:
2804  {
2805  switch (shape)
2806  {
2807  case LibUtilities::eSegment:
2808  {
2809  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
2810  LibUtilities::BasisKey bkey(LibUtilities::eOrtho_A, nummodes, pkey);
2811 
2812  returnval.push_back(bkey);
2813  }
2814  break;
2815  case LibUtilities::eTriangle:
2816  {
2817  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
2818  LibUtilities::BasisKey bkey(LibUtilities::eOrtho_A, nummodes, pkey);
2819 
2820  returnval.push_back(bkey);
2821 
2822  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussRadauMAlpha1Beta0);
2823  LibUtilities::BasisKey bkey1(LibUtilities::eOrtho_B, nummodes, pkey1);
2824 
2825  returnval.push_back(bkey1);
2826  }
2827  break;
2828  case LibUtilities::eQuadrilateral:
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  returnval.push_back(bkey);
2835  }
2836  break;
2837  case LibUtilities::eTetrahedron:
2838  {
2839  const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussLobattoLegendre);
2840  LibUtilities::BasisKey bkey(LibUtilities::eOrtho_A, nummodes, pkey);
2841 
2842  returnval.push_back(bkey);
2843 
2844  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussRadauMAlpha1Beta0);
2845  LibUtilities::BasisKey bkey1(LibUtilities::eOrtho_B, nummodes, pkey1);
2846 
2847  returnval.push_back(bkey1);
2848 
2849  const LibUtilities::PointsKey pkey2(nummodes, LibUtilities::eGaussRadauMAlpha2Beta0);
2850  LibUtilities::BasisKey bkey2(LibUtilities::eOrtho_C, nummodes, pkey2);
2851  }
2852  break;
2853  default:
2854  {
2855  ASSERTL0(false,"Expansion not defined in switch for this shape");
2856  }
2857  break;
2858  }
2859  }
2860  break;
2861 
2862  case eGLL_Lagrange_SEM:
2863  {
2864  switch (shape)
2865  {
2866  case LibUtilities::eSegment:
2867  {
2868  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussLobattoLegendre);
2869  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
2870 
2871  returnval.push_back(bkey);
2872  }
2873  break;
2874  case LibUtilities::eQuadrilateral:
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  }
2882  break;
2883  case LibUtilities::eHexahedron:
2884  {
2885  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussLobattoLegendre);
2886  LibUtilities::BasisKey bkey(LibUtilities::eGLL_Lagrange, nummodes, pkey);
2887 
2888  returnval.push_back(bkey);
2889  returnval.push_back(bkey);
2890  returnval.push_back(bkey);
2891  }
2892  break;
2893  default:
2894  {
2895  ASSERTL0(false,"Expansion not defined in switch for this shape");
2896  }
2897  break;
2898  }
2899  }
2900  break;
2901 
2902 
2903  case eFourier:
2904  {
2905  switch (shape)
2906  {
2907  case LibUtilities::eSegment:
2908  {
2909  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
2910  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
2911  returnval.push_back(bkey);
2912  }
2913  break;
2914  case LibUtilities::eQuadrilateral:
2915  {
2916  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
2917  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
2918  returnval.push_back(bkey);
2919  returnval.push_back(bkey);
2920  }
2921  break;
2922  case LibUtilities::eHexahedron:
2923  {
2924  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
2925  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
2926  returnval.push_back(bkey);
2927  returnval.push_back(bkey);
2928  returnval.push_back(bkey);
2929  }
2930  break;
2931  default:
2932  {
2933  ASSERTL0(false,"Expansion not defined in switch for this shape");
2934  }
2935  break;
2936  }
2937  }
2938  break;
2939 
2940 
2941  case eFourierSingleMode:
2942  {
2943  switch (shape)
2944  {
2945  case LibUtilities::eSegment:
2946  {
2947  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
2948  LibUtilities::BasisKey bkey(LibUtilities::eFourierSingleMode, nummodes, pkey);
2949  returnval.push_back(bkey);
2950  }
2951  break;
2952  case LibUtilities::eQuadrilateral:
2953  {
2954  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
2955  LibUtilities::BasisKey bkey(LibUtilities::eFourierSingleMode, nummodes, pkey);
2956  returnval.push_back(bkey);
2957  returnval.push_back(bkey);
2958  }
2959  break;
2960  case LibUtilities::eHexahedron:
2961  {
2962  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
2963  LibUtilities::BasisKey bkey(LibUtilities::eFourierSingleMode, nummodes, pkey);
2964  returnval.push_back(bkey);
2965  returnval.push_back(bkey);
2966  returnval.push_back(bkey);
2967  }
2968  break;
2969  default:
2970  {
2971  ASSERTL0(false,"Expansion not defined in switch for this shape");
2972  }
2973  break;
2974  }
2975  }
2976  break;
2977 
2978  case eFourierHalfModeRe:
2979  {
2980  switch (shape)
2981  {
2982  case LibUtilities::eSegment:
2983  {
2984  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
2985  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeRe, nummodes, pkey);
2986  returnval.push_back(bkey);
2987  }
2988  break;
2989  case LibUtilities::eQuadrilateral:
2990  {
2991  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
2992  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeRe, nummodes, pkey);
2993  returnval.push_back(bkey);
2994  returnval.push_back(bkey);
2995  }
2996  break;
2997  case LibUtilities::eHexahedron:
2998  {
2999  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3000  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeRe, nummodes, pkey);
3001  returnval.push_back(bkey);
3002  returnval.push_back(bkey);
3003  returnval.push_back(bkey);
3004  }
3005  break;
3006  default:
3007  {
3008  ASSERTL0(false,"Expansion not defined in switch for this shape");
3009  }
3010  break;
3011  }
3012  }
3013  break;
3014 
3015  case eFourierHalfModeIm:
3016  {
3017  switch (shape)
3018  {
3019  case LibUtilities::eSegment:
3020  {
3021  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3022  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeIm, nummodes, pkey);
3023  returnval.push_back(bkey);
3024  }
3025  break;
3026  case LibUtilities::eQuadrilateral:
3027  {
3028  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3029  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeIm, nummodes, pkey);
3030  returnval.push_back(bkey);
3031  returnval.push_back(bkey);
3032  }
3033  break;
3034  case LibUtilities::eHexahedron:
3035  {
3036  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierSingleModeSpaced);
3037  LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeIm, nummodes, pkey);
3038  returnval.push_back(bkey);
3039  returnval.push_back(bkey);
3040  returnval.push_back(bkey);
3041  }
3042  break;
3043  default:
3044  {
3045  ASSERTL0(false,"Expansion not defined in switch for this shape");
3046  }
3047  break;
3048  }
3049  }
3050  break;
3051 
3052  case eChebyshev:
3053  {
3054  switch (shape)
3055  {
3056  case LibUtilities::eSegment:
3057  {
3058  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussChebyshev);
3059  LibUtilities::BasisKey bkey(LibUtilities::eChebyshev, nummodes, pkey);
3060  returnval.push_back(bkey);
3061  }
3062  break;
3063  case LibUtilities::eQuadrilateral:
3064  {
3065  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussChebyshev);
3066  LibUtilities::BasisKey bkey(LibUtilities::eChebyshev, nummodes, pkey);
3067  returnval.push_back(bkey);
3068  returnval.push_back(bkey);
3069  }
3070  break;
3071  case LibUtilities::eHexahedron:
3072  {
3073  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussChebyshev);
3074  LibUtilities::BasisKey bkey(LibUtilities::eChebyshev, nummodes, pkey);
3075  returnval.push_back(bkey);
3076  returnval.push_back(bkey);
3077  returnval.push_back(bkey);
3078  }
3079  break;
3080  default:
3081  {
3082  ASSERTL0(false,"Expansion not defined in switch for this shape");
3083  }
3084  break;
3085  }
3086  }
3087  break;
3088 
3089  case eFourierChebyshev:
3090  {
3091  switch (shape)
3092  {
3093  case LibUtilities::eQuadrilateral:
3094  {
3095  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3096  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3097  returnval.push_back(bkey);
3098 
3099  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussGaussChebyshev);
3100  LibUtilities::BasisKey bkey1(LibUtilities::eChebyshev, nummodes, pkey1);
3101  returnval.push_back(bkey1);
3102  }
3103  break;
3104  default:
3105  {
3106  ASSERTL0(false,"Expansion not defined in switch for this shape");
3107  }
3108  break;
3109  }
3110  }
3111  break;
3112 
3113  case eChebyshevFourier:
3114  {
3115  switch (shape)
3116  {
3117  case LibUtilities::eQuadrilateral:
3118  {
3119  const LibUtilities::PointsKey pkey1(nummodes, LibUtilities::eGaussGaussChebyshev);
3120  LibUtilities::BasisKey bkey1(LibUtilities::eChebyshev, nummodes, pkey1);
3121  returnval.push_back(bkey1);
3122 
3123  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3124  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3125  returnval.push_back(bkey);
3126  }
3127  break;
3128  default:
3129  {
3130  ASSERTL0(false,"Expansion not defined in switch for this shape");
3131  }
3132  break;
3133  }
3134  }
3135  break;
3136 
3137  case eFourierModified:
3138  {
3139  switch (shape)
3140  {
3141  case LibUtilities::eQuadrilateral:
3142  {
3143  const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eFourierEvenlySpaced);
3144  LibUtilities::BasisKey bkey(LibUtilities::eFourier, nummodes, pkey);
3145  returnval.push_back(bkey);
3146 
3147  const LibUtilities::PointsKey pkey1(nummodes+1, LibUtilities::eGaussLobattoLegendre);
3148  LibUtilities::BasisKey bkey1(LibUtilities::eModified_A, nummodes, pkey1);
3149  returnval.push_back(bkey1);
3150  }
3151  break;
3152  default:
3153  {
3154  ASSERTL0(false,"Expansion not defined in switch for this shape");
3155  }
3156  break;
3157  }
3158  }
3159  break;
3160 
3161  default:
3162  {
3163  ASSERTL0(false,"Expansion type not defined");
3164  }
3165  break;
3166 
3167  }
3168 
3169  return returnval;
3170  }
3171 
3172  /**
3173  *
3174  */
3175  LibUtilities::BasisKeyVector
3176  MeshGraph::DefineBasisKeyFromExpansionTypeHomo(
3177  GeometrySharedPtr in,
3178  ExpansionType type_x,
3179  ExpansionType type_y,
3180  ExpansionType type_z,
3181  const int nummodes_x,
3182  const int nummodes_y,
3183  const int nummodes_z)
3184  {
3185  LibUtilities::BasisKeyVector returnval;
3186 
3187  LibUtilities::ShapeType shape = in->GetShapeType();
3188 
3189  switch (shape)
3190  {
3191  case LibUtilities::eSegment:
3192  {
3193  ASSERTL0(false,"Homogeneous expansion not defined for this shape");
3194  }
3195  break;
3196 
3197  case LibUtilities::eQuadrilateral:
3198  {
3199  ASSERTL0(false,"Homogeneous expansion not defined for this shape");
3200  }
3201  break;
3202 
3203  case LibUtilities::eHexahedron:
3204  {
3205  switch(type_x)
3206  {
3207  case eFourier:
3208  {
3209  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eFourierEvenlySpaced);
3210  LibUtilities::BasisKey bkey1(LibUtilities::eFourier,nummodes_x,pkey1);
3211  returnval.push_back(bkey1);
3212  }
3213  break;
3214 
3215  case eFourierSingleMode:
3216  {
3217  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eFourierSingleModeSpaced);
3218  LibUtilities::BasisKey bkey1(LibUtilities::eFourierSingleMode,nummodes_x,pkey1);
3219  returnval.push_back(bkey1);
3220  }
3221  break;
3222 
3223  case eFourierHalfModeRe:
3224  {
3225  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eFourierSingleModeSpaced);
3226  LibUtilities::BasisKey bkey1(LibUtilities::eFourierHalfModeRe,nummodes_x,pkey1);
3227  returnval.push_back(bkey1);
3228  }
3229  break;
3230 
3231  case eFourierHalfModeIm:
3232  {
3233  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eFourierSingleModeSpaced);
3234  LibUtilities::BasisKey bkey1(LibUtilities::eFourierHalfModeIm,nummodes_x,pkey1);
3235  returnval.push_back(bkey1);
3236  }
3237  break;
3238 
3239 
3240  case eChebyshev:
3241  {
3242  const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eGaussGaussChebyshev);
3243  LibUtilities::BasisKey bkey1(LibUtilities::eChebyshev,nummodes_x,pkey1);
3244  returnval.push_back(bkey1);
3245  }
3246  break;
3247 
3248 
3249 
3250  default:
3251  {
3252  ASSERTL0(false,"Homogeneous expansion can be of Fourier or Chebyshev type only");
3253  }
3254  break;
3255  }
3256 
3257 
3258  switch(type_y)
3259  {
3260  case eFourier:
3261  {
3262  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eFourierEvenlySpaced);
3263  LibUtilities::BasisKey bkey2(LibUtilities::eFourier,nummodes_y,pkey2);
3264  returnval.push_back(bkey2);
3265  }
3266  break;
3267 
3268 
3269  case eFourierSingleMode:
3270  {
3271  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eFourierSingleModeSpaced);
3272  LibUtilities::BasisKey bkey2(LibUtilities::eFourierSingleMode,nummodes_y,pkey2);
3273  returnval.push_back(bkey2);
3274  }
3275  break;
3276 
3277  case eFourierHalfModeRe:
3278  {
3279  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eFourierSingleModeSpaced);
3280  LibUtilities::BasisKey bkey2(LibUtilities::eFourierHalfModeRe,nummodes_y,pkey2);
3281  returnval.push_back(bkey2);
3282  }
3283  break;
3284 
3285  case eFourierHalfModeIm:
3286  {
3287  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eFourierSingleModeSpaced);
3288  LibUtilities::BasisKey bkey2(LibUtilities::eFourierHalfModeIm,nummodes_y,pkey2);
3289  returnval.push_back(bkey2);
3290  }
3291  break;
3292 
3293  case eChebyshev:
3294  {
3295  const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eGaussGaussChebyshev);
3296  LibUtilities::BasisKey bkey2(LibUtilities::eChebyshev,nummodes_y,pkey2);
3297  returnval.push_back(bkey2);
3298  }
3299  break;
3300 
3301  default:
3302  {
3303  ASSERTL0(false,"Homogeneous expansion can be of Fourier or Chebyshev type only");
3304  }
3305  break;
3306  }
3307 
3308  switch(type_z)
3309  {
3310  case eFourier:
3311  {
3312  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eFourierEvenlySpaced);
3313  LibUtilities::BasisKey bkey3(LibUtilities::eFourier,nummodes_z,pkey3);
3314  returnval.push_back(bkey3);
3315  }
3316  break;
3317 
3318  case eFourierSingleMode:
3319  {
3320  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eFourierSingleModeSpaced);
3321  LibUtilities::BasisKey bkey3(LibUtilities::eFourierSingleMode,nummodes_z,pkey3);
3322  returnval.push_back(bkey3);
3323  }
3324  break;
3325 
3326  case eFourierHalfModeRe:
3327  {
3328  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eFourierSingleModeSpaced);
3329  LibUtilities::BasisKey bkey3(LibUtilities::eFourierHalfModeRe,nummodes_z,pkey3);
3330  returnval.push_back(bkey3);
3331  }
3332  break;
3333 
3334  case eFourierHalfModeIm:
3335  {
3336  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eFourierSingleModeSpaced);
3337  LibUtilities::BasisKey bkey3(LibUtilities::eFourierHalfModeIm,nummodes_z,pkey3);
3338  returnval.push_back(bkey3);
3339  }
3340  break;
3341 
3342  case eChebyshev:
3343  {
3344  const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eGaussGaussChebyshev);
3345  LibUtilities::BasisKey bkey3(LibUtilities::eChebyshev,nummodes_z,pkey3);
3346  returnval.push_back(bkey3);
3347  }
3348  break;
3349 
3350  default:
3351  {
3352  ASSERTL0(false,"Homogeneous expansion can be of Fourier or Chebyshev type only");
3353  }
3354  break;
3355  }
3356  }
3357  break;
3358 
3359  case LibUtilities::eTriangle:
3360  {
3361  ASSERTL0(false,"Homogeneous expansion not defined for this shape");
3362  }
3363  break;
3364 
3365  case LibUtilities::eTetrahedron:
3366  {
3367  ASSERTL0(false,"Homogeneous expansion not defined for this shape");
3368  }
3369  break;
3370 
3371  default:
3372  ASSERTL0(false,"Expansion not defined in switch for this shape");
3373  break;
3374  }
3375 
3376  return returnval;
3377  }
3378 
3379 
3380  /**
3381  *
3382  */
3383  PointGeomSharedPtr MeshGraph::AddVertex(NekDouble x, NekDouble y, NekDouble z)
3384  {
3385  unsigned int nextId = m_vertSet.rbegin()->first + 1;
3386  PointGeomSharedPtr vert(MemoryManager<PointGeom>::AllocateSharedPtr(m_spaceDimension, nextId, x, y, z));
3387  m_vertSet[nextId] = vert;
3388  return vert;
3389  }
3390 
3391 
3392  /**
3393  *
3394  */
3395  SegGeomSharedPtr MeshGraph::AddEdge(PointGeomSharedPtr v0, PointGeomSharedPtr v1,
3396  CurveSharedPtr curveDefinition)
3397  {
3398  PointGeomSharedPtr vertices[] = {v0, v1};
3399  SegGeomSharedPtr edge;
3400  int edgeId = m_segGeoms.rbegin()->first + 1;
3401 
3402  if( curveDefinition )
3403  {
3404  edge = MemoryManager<SegGeom>::AllocateSharedPtr(edgeId, m_spaceDimension, vertices, curveDefinition);
3405  }
3406  else
3407  {
3408  edge = MemoryManager<SegGeom>::AllocateSharedPtr(edgeId, m_spaceDimension, vertices);
3409  }
3410  m_segGeoms[edgeId] = edge;
3411  return edge;
3412  }
3413 
3414 
3415  /**
3416  *
3417  */
3418  TriGeomSharedPtr MeshGraph::AddTriangle(SegGeomSharedPtr edges[], StdRegions::Orientation orient[])
3419  {
3420  int indx = m_triGeoms.rbegin()->first + 1;
3421  TriGeomSharedPtr trigeom(MemoryManager<TriGeom>::AllocateSharedPtr(indx, edges, orient));
3422  trigeom->SetGlobalID(indx);
3423 
3424  m_triGeoms[indx] = trigeom;
3425 
3426  return trigeom;
3427  }
3428 
3429 
3430  /**
3431  *
3432  */
3433  QuadGeomSharedPtr MeshGraph::AddQuadrilateral(SegGeomSharedPtr edges[], StdRegions::Orientation orient[])
3434  {
3435  int indx = m_quadGeoms.rbegin()->first + 1;
3436  QuadGeomSharedPtr quadgeom(MemoryManager<QuadGeom>::AllocateSharedPtr(indx, edges, orient));
3437  quadgeom->SetGlobalID(indx);
3438 
3439  m_quadGeoms[indx] = quadgeom;
3440  return quadgeom;
3441  }
3442 
3443 
3444  /**
3445  *
3446  */
3447  PrismGeomSharedPtr MeshGraph::AddPrism(TriGeomSharedPtr tfaces[PrismGeom::kNtfaces],
3448  QuadGeomSharedPtr qfaces[PrismGeom::kNqfaces])
3449  {
3450  // Setting the orientation is disabled in the reader. Why?
3451  Geometry2DSharedPtr faces[] = { qfaces[0], tfaces[0], qfaces[1], tfaces[1], qfaces[2] };
3452  unsigned int index = m_prismGeoms.rbegin()->first + 1;
3453  PrismGeomSharedPtr prismgeom(MemoryManager<PrismGeom>::AllocateSharedPtr(faces));
3454  prismgeom->SetGlobalID(index);
3455 
3456  m_prismGeoms[index] = prismgeom;
3457  return prismgeom;
3458  }
3459 
3460 
3461  /**
3462  *
3463  */
3464  TetGeomSharedPtr MeshGraph::AddTetrahedron(TriGeomSharedPtr tfaces[TetGeom::kNtfaces])
3465  {
3466  unsigned int index = m_tetGeoms.rbegin()->first + 1;
3467  TetGeomSharedPtr tetgeom(MemoryManager<TetGeom>::AllocateSharedPtr(tfaces));
3468  tetgeom->SetGlobalID(index);
3469 
3470  m_tetGeoms[index] = tetgeom;
3471  return tetgeom;
3472  }
3473 
3474 
3475  /**
3476  *
3477  */
3478  PyrGeomSharedPtr MeshGraph::AddPyramid(TriGeomSharedPtr tfaces[PyrGeom::kNtfaces],
3479  QuadGeomSharedPtr qfaces[PyrGeom::kNqfaces])
3480  {
3481  Geometry2DSharedPtr faces[] = { qfaces[0], tfaces[0], tfaces[1], tfaces[2], tfaces[3] };
3482  unsigned int index = m_pyrGeoms.rbegin()->first + 1;
3483 
3484  PyrGeomSharedPtr pyrgeom(MemoryManager<PyrGeom>::AllocateSharedPtr(faces));
3485  pyrgeom->SetGlobalID(index);
3486 
3487  m_pyrGeoms[index] = pyrgeom;
3488  return pyrgeom;
3489  }
3490 
3491 
3492  /**
3493  *
3494  */
3495  HexGeomSharedPtr MeshGraph::AddHexahedron(QuadGeomSharedPtr qfaces[HexGeom::kNqfaces])
3496  {
3497  unsigned int index = m_hexGeoms.rbegin()->first + 1;
3498  HexGeomSharedPtr hexgeom(MemoryManager<HexGeom>::AllocateSharedPtr(qfaces));
3499  hexgeom->SetGlobalID(index);
3500  m_hexGeoms[index] = hexgeom;
3501  return hexgeom;
3502  }
3503 
3504 
3505  /**
3506  * Generate a single vector of Expansion structs mapping global element
3507  * ID to a corresponding Geometry shared pointer and basis key.
3508  *
3509  * Expansion map ensures elements which appear in multiple composites
3510  * within the domain are only listed once.
3511  */
3512  ExpansionMapShPtr MeshGraph::SetUpExpansionMap(void)
3513  {
3514  ExpansionMapShPtr returnval;
3515  returnval = MemoryManager<ExpansionMap>::AllocateSharedPtr();
3516 
3517  for(int d = 0; d < m_domain.size(); ++d)
3518  {
3519  CompositeMap::const_iterator compIter;
3520 
3521  for (compIter = m_domain[d].begin(); compIter != m_domain[d].end(); ++compIter)
3522  {
3523  GeometryVector::const_iterator x;
3524  for (x = compIter->second->begin(); x != compIter->second->end(); ++x)
3525  {
3526  LibUtilities::BasisKeyVector def;
3527  ExpansionShPtr expansionElementShPtr =
3528  MemoryManager<Expansion>::AllocateSharedPtr(*x, def);
3529  int id = (*x)->GetGlobalID();
3530  (*returnval)[id] = expansionElementShPtr;
3531  }
3532  }
3533  }
3534 
3535  return returnval;
3536  }
3537  }; //end of namespace
3538 }; //end of namespace
Generated on Wed Mar 4 2015 16:33:15 for Nektar++ by   doxygen 1.8.1.2