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