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