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