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