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