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