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