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