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