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