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