#include <InputNek.h>

Inheritance diagram for Nektar::Utilities::InputNek:

Public Member Functions
	InputNek (NekMeshUtils::MeshSharedPtr p_m)

virtual	~InputNek ()

virtual void	Process ()
	Processes Nektar file format. More...

Public Member Functions inherited from Nektar::NekMeshUtils::InputModule
NEKMESHUTILS_EXPORT	InputModule (MeshSharedPtr p_m)

NEKMESHUTILS_EXPORT void	OpenStream ()
	Open a file for input. More...

Public Member Functions inherited from Nektar::NekMeshUtils::Module
NEKMESHUTILS_EXPORT	Module (MeshSharedPtr p_m)

NEKMESHUTILS_EXPORT void	RegisterConfig (std::string key, std::string value=std::string())
	Register a configuration option with a module. More...

NEKMESHUTILS_EXPORT void	PrintConfig ()
	Print out all configuration options for a module. More...

NEKMESHUTILS_EXPORT void	SetDefaults ()
	Sets default configuration options for those which have not been set. More...

NEKMESHUTILS_EXPORT MeshSharedPtr	GetMesh ()

virtual NEKMESHUTILS_EXPORT void	ProcessVertices ()
	Extract element vertices. More...

virtual NEKMESHUTILS_EXPORT void	ProcessEdges (bool ReprocessEdges=true)
	Extract element edges. More...

virtual NEKMESHUTILS_EXPORT void	ProcessFaces (bool ReprocessFaces=true)
	Extract element faces. More...

virtual NEKMESHUTILS_EXPORT void	ProcessElements ()
	Generate element IDs. More...

virtual NEKMESHUTILS_EXPORT void	ProcessComposites ()
	Generate composites. More...

virtual NEKMESHUTILS_EXPORT void	ClearElementLinks ()

Static Public Member Functions
static NekMeshUtils::ModuleSharedPtr	create (NekMeshUtils::MeshSharedPtr m)
	Creates an instance of this class. More...

Static Public Attributes
static NekMeshUtils::ModuleKey	className
	ModuleKey for class. More...

Private Member Functions
void	LoadHOSurfaces ()

int	GetNnodes (LibUtilities::ShapeType elType)

Private Attributes
std::map< std::string, std::pair< NekCurve, std::string > >	curveTags

std::map< std::string, NekMeshUtils::HOSurfSet >	hoData

std::map< int, int >	hoMap

Additional Inherited Members
Protected Member Functions inherited from Nektar::NekMeshUtils::InputModule
NEKMESHUTILS_EXPORT void	PrintSummary ()
	Print summary of elements. More...

Protected Member Functions inherited from Nektar::NekMeshUtils::Module
NEKMESHUTILS_EXPORT void	ReorderPrisms (PerMap &perFaces)
	Reorder node IDs so that prisms and tetrahedra are aligned correctly. More...

NEKMESHUTILS_EXPORT void	PrismLines (int prism, PerMap &perFaces, std::set< int > &prismsDone, std::vector< ElementSharedPtr > &line)

Protected Attributes inherited from Nektar::NekMeshUtils::InputModule
io::filtering_istream	m_mshFile
	Input stream. More...

std::ifstream	m_mshFileStream
	Input stream. More...

Protected Attributes inherited from Nektar::NekMeshUtils::Module
MeshSharedPtr	m_mesh
	Mesh object. More...

std::map< std::string, ConfigOption >	m_config
	List of configuration values. More...

Detailed Description

Converter class for Nektar session files.

Definition at line 59 of file InputNek.h.

Constructor & Destructor Documentation

◆ InputNek()

Nektar::Utilities::InputNek::InputNek ( NekMeshUtils::MeshSharedPtr p_m )

Definition at line 59 of file InputNek.cpp.

References Nektar::NekMeshUtils::Module::m_config.

                                   : InputModule(m)
 {
     m_config["scalar"] = ConfigOption(
         true, "0", "If defined then assume input rea is for scalar problem");
 }

◆ ~InputNek()

Nektar::Utilities::InputNek::~InputNek ( )

virtual

Definition at line 65 of file InputNek.cpp.

66 {

67 }

Member Function Documentation

◆ create()

static NekMeshUtils::ModuleSharedPtr Nektar::Utilities::InputNek::create ( NekMeshUtils::MeshSharedPtr m )

inlinestatic

Creates an instance of this class.

Definition at line 67 of file InputNek.h.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

     {
         return MemoryManager<InputNek>::AllocateSharedPtr(m);
     }

◆ GetNnodes()

int Nektar::Utilities::InputNek::GetNnodes ( LibUtilities::ShapeType InputNekEntity )

private

This routine aids the reading of Nektar files only; it returns the number of nodes for a given entity typw.

Definition at line 1114 of file InputNek.cpp.

References Nektar::LibUtilities::eHexahedron, Nektar::LibUtilities::ePoint, Nektar::LibUtilities::ePrism, Nektar::LibUtilities::ePyramid, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eSegment, Nektar::LibUtilities::eTetrahedron, and Nektar::LibUtilities::eTriangle.

Referenced by Process().

 {
     int nNodes = 0;
 
     switch (InputNekEntity)
     {
         case LibUtilities::ePoint:
             nNodes = 1;
             break;
         case LibUtilities::eSegment:
             nNodes = 2;
             break;
         case LibUtilities::eTriangle:
             nNodes = 3;
             break;
         case LibUtilities::eQuadrilateral:
             nNodes = 4;
             break;
         case LibUtilities::eTetrahedron:
             nNodes = 4;
             break;
         case LibUtilities::ePyramid:
             nNodes = 5;
             break;
         case LibUtilities::ePrism:
             nNodes = 6;
             break;
         case LibUtilities::eHexahedron:
             nNodes = 8;
             break;
         default:
             cerr << "unknown Nektar element type" << endl;
     }
 
     return nNodes;
 }

◆ LoadHOSurfaces()

void Nektar::Utilities::InputNek::LoadHOSurfaces ( )

private

Load high order surface information from hsf file.

Definition at line 950 of file InputNek.cpp.

References curveTags, Nektar::Utilities::eFile, Nektar::LibUtilities::eNodalTriElec, hoData, hoMap, Nektar::NekMeshUtils::Module::m_mesh, class_topology::Node, and Nektar::LibUtilities::PointsManager().

Referenced by Process().

 {
     int nodeId = m_mesh->GetNumEntities();
 
     for (auto &it : curveTags)
     {
         ifstream hsf;
         string line, fileName = it.second.second;
         size_t pos;
         int N, Nface, dot;
 
         if (it.second.first != eFile)
         {
             continue;
         }
 
         // Replace fro extension with hsf.
         dot      = fileName.find_last_of('.');
         fileName = fileName.substr(0, dot);
         fileName += ".hsf";
 
         // Open hsf file.
         hsf.open(fileName.c_str());
         if (!hsf.is_open())
         {
             cerr << "Could not open surface file " << fileName << endl;
             abort();
         }
 
         // Read in header line; determine element order, number of faces
         // from this line.
         getline(hsf, line);
         pos = line.find("=");
         if (pos == string::npos)
         {
             cerr << "hsf header error: cannot read number of "
                  << "nodal points." << endl;
             abort();
         }
         line = line.substr(pos + 1);
         stringstream ss(line);
         ss >> N;
 
         pos = line.find("=");
         if (pos == string::npos)
         {
             cerr << "hsf header error: cannot read number of "
                  << "faces." << endl;
             abort();
         }
         line = line.substr(pos + 1);
         ss.clear();
         ss.str(line);
         ss >> Nface;
 
         int Ntot = N * (N + 1) / 2;
 
         // Skip a line, then read in r,s positions inside the next
         // comments.
         Array<OneD, NekDouble> r(Ntot), s(Ntot);
         getline(hsf, line);
 
         for (int i = 0; i < 2; ++i)
         {
             string word;
 
             getline(hsf, line);
             ss.clear();
             ss.str(line);
             ss >> word;
 
             if (word != "#")
             {
                 cerr << "hsf header error: cannot read in "
                      << "r/s points" << endl;
                 abort();
             }
 
             for (int j = 0; j < Ntot; ++j)
             {
                 ss >> (i == 0 ? r[j] : s[j]);
             }
         }
 
         // Generate electrostatic points so that re-mapping array can
         // be constructed.
         Array<OneD, NekDouble> rp(Ntot), sp(Ntot);
         LibUtilities::PointsKey elec(N, LibUtilities::eNodalTriElec);
         LibUtilities::PointsManager()[elec]->GetPoints(rp, sp);
 
         // Expensively construct remapping array nodemap. This will
         // map nodal ordering from hsf order to Nektar++ ordering
         // (i.e. vertices followed by edges followed by interior
         // points.)
         for (int i = 0; i < Ntot; ++i)
         {
             for (int j = 0; j < Ntot; ++j)
             {
                 if (fabs(r[i] - rp[j]) < 1e-5 && fabs(s[i] - sp[j]) < 1e-5)
                 {
                     hoMap[i] = j;
                     break;
                 }
             }
         }
 
         // Skip variables line
         getline(hsf, line);
 
         // Read in nodal points for each face.
         map<int, vector<NodeSharedPtr> > faceMap;
         for (int i = 0; i < Nface; ++i)
         {
             getline(hsf, line);
             vector<NodeSharedPtr> faceNodes(Ntot);
             for (int j = 0; j < Ntot; ++j, ++nodeId)
             {
                 double x, y, z;
                 getline(hsf, line);
                 ss.clear();
                 ss.str(line);
                 ss >> x >> y >> z;
                 faceNodes[j] = NodeSharedPtr(new Node(nodeId, x, y, z));
             }
             // Skip over tecplot connectivity information.
             for (int j = 0; j < (N - 1) * (N - 1); ++j)
             {
                 getline(hsf, line);
             }
             faceMap[i] = faceNodes;
         }
 
         // Finally, read in connectivity information to set up after
         // reading rea file.
         getline(hsf, line);
         for (int i = 0; i < Nface; ++i)
         {
             string tmp;
             int fid;
             vector<int> nodeIds(3);
 
             getline(hsf, line);
             ss.clear();
             ss.str(line);
             ss >> tmp >> fid >> nodeIds[0] >> nodeIds[1] >> nodeIds[2];
 
             if (tmp != "#")
             {
                 cerr << "Unable to read hsf connectivity information." << endl;
                 abort();
             }
 
             hoData[it.first].insert(
                 HOSurfSharedPtr(new HOSurf(nodeIds, faceMap[i])));
         }
 
         hsf.close();
     }
 }

◆ Process()

void Nektar::Utilities::InputNek::Process ( )

virtual

Processes Nektar file format.

Nektar sessions are defined by rea files, and contain sections defining a DNS simulation in a specific order. The converter only reads mesh information, curve information if it exists and boundary information.

Parameters

pFilename Filename of Nektar session file to read.

Implements Nektar::NekMeshUtils::Module.

Definition at line 78 of file InputNek.cpp.

 {
     // Open the file stream.
     OpenStream();
 
     string line, word;
     int nParam, nElements, nCurves;
     int i, j, k, nodeCounter = 0;
     int nComposite           = 0;
     LibUtilities::ShapeType elType;
     double vertex[3][8];
     map<LibUtilities::ShapeType, int> domainComposite;
     map<LibUtilities::ShapeType, vector<vector<NodeSharedPtr> > > elNodes;
     map<LibUtilities::ShapeType, vector<int> > elIds;
     std::unordered_map<int, int> elMap;
     vector<LibUtilities::ShapeType> elmOrder;
 
     bool scalar = m_config["scalar"].as<bool>();
 
     // Set up vector of processing orders.
     elmOrder.push_back(LibUtilities::eSegment);
     elmOrder.push_back(LibUtilities::eTriangle);
     elmOrder.push_back(LibUtilities::eQuadrilateral);
     elmOrder.push_back(LibUtilities::ePrism);
     elmOrder.push_back(LibUtilities::ePyramid);
     elmOrder.push_back(LibUtilities::eTetrahedron);
     elmOrder.push_back(LibUtilities::eHexahedron);
 
     m_mesh->m_expDim   = 0;
     m_mesh->m_spaceDim = 0;
 
     if (m_mesh->m_verbose)
     {
         cout << "InputNek: Start reading file..." << endl;
     }
 
     // -- Read in parameters.
 
     // Ignore first 3 lines. 4th line contains number of parameters.
     for (i = 0; i < 4; ++i)
     {
         getline(m_mshFile, line);
     }
 
     stringstream s(line);
     s >> nParam;
 
     for (i = 0; i < nParam; ++i)
     {
         string tmp1, tmp2;
         getline(m_mshFile, line);
         s.str(line);
         s >> tmp1 >> tmp2;
     }
 
     // -- Read in passive scalars (ignore)
     getline(m_mshFile, line);
     s.clear();
     s.str(line);
     s >> j;
     for (i = 0; i < j; ++i)
     {
         getline(m_mshFile, line);
     }
 
     // -- Read in logical switches (ignore)
     getline(m_mshFile, line);
     s.clear();
     s.str(line);
     s >> j;
     for (i = 0; i < j; ++i)
     {
         getline(m_mshFile, line);
     }
 
     // -- Read in mesh data.
 
     // First hunt for MESH tag
     bool foundMesh = false;
     while (!m_mshFile.eof())
     {
         getline(m_mshFile, line);
         if (line.find("MESH") != string::npos)
         {
             foundMesh = true;
             break;
         }
     }
 
     if (!foundMesh)
     {
         cerr << "Couldn't find MESH tag inside file." << endl;
         abort();
     }
 
     // Now read in number of elements and space dimension.
     getline(m_mshFile, line);
     s.clear();
     s.str(line);
     s >> nElements >> m_mesh->m_expDim;
     m_mesh->m_spaceDim = m_mesh->m_expDim;
 
     // Set up field names.
     m_mesh->m_fields.push_back("u");
 
     if (!scalar)
     {
         m_mesh->m_fields.push_back("v");
         if (m_mesh->m_spaceDim > 2)
         {
             m_mesh->m_fields.push_back("w");
         }
         m_mesh->m_fields.push_back("p");
     }
 
     // Loop over and create elements.
     for (i = 0; i < nElements; ++i)
     {
         getline(m_mshFile, line);
 
         if (m_mesh->m_expDim == 2)
         {
             if (line.find("Qua") != string::npos ||
                 line.find("qua") != string::npos)
             {
                 elType = LibUtilities::eQuadrilateral;
             }
             else
             {
                 // Default element type in 2D is triangle.
                 elType = LibUtilities::eTriangle;
             }
         }
         else
         {
             if (line.find("Tet") != string::npos ||
                 line.find("tet") != string::npos)
             {
                 elType = LibUtilities::eTetrahedron;
             }
             else if (line.find("Hex") != string::npos ||
                      line.find("hex") != string::npos)
             {
                 elType = LibUtilities::eHexahedron;
             }
             else if (line.find("Prism") != string::npos ||
                      line.find("prism") != string::npos)
             {
                 elType = LibUtilities::ePrism;
             }
             else if (line.find("Pyr") != string::npos ||
                      line.find("pyr") != string::npos)
             {
                 elType = LibUtilities::ePyramid;
             }
             else if (line.find("Qua") != string::npos ||
                      line.find("qua") != string::npos)
             {
                 elType = LibUtilities::eQuadrilateral;
             }
             else
             {
                 // Default element type in 2D is tetrahedron.
                 elType = LibUtilities::eTetrahedron;
             }
         }
 
         // Read in number of vertices for element type.
         const int nNodes = GetNnodes(elType);
 
         for (j = 0; j < m_mesh->m_expDim; ++j)
         {
             getline(m_mshFile, line);
             s.clear();
             s.str(line);
             for (k = 0; k < nNodes; ++k)
             {
                 s >> vertex[j][k];
             }
         }
 
         // Zero co-ordinates bigger than expansion dimension.
         for (j = m_mesh->m_expDim; j < 3; ++j)
         {
             for (k = 0; k < nNodes; ++k)
             {
                 vertex[j][k] = 0.0;
             }
         }
 
         // Nektar meshes do not contain a unique list of nodes, so this
         // block constructs a unique set so that elements can be created
         // with unique nodes.
         vector<NodeSharedPtr> nodeList;
         for (k = 0; k < nNodes; ++k)
         {
             NodeSharedPtr n = std::shared_ptr<Node>(new Node(
                 nodeCounter++, vertex[0][k], vertex[1][k], vertex[2][k]));
             nodeList.push_back(n);
         }
 
         elNodes[elType].push_back(nodeList);
         elIds[elType].push_back(i);
     }
 
     int reorderedId = 0;
     nodeCounter     = 0;
 
     for (i = 0; i < elmOrder.size(); ++i)
     {
         LibUtilities::ShapeType elType      = elmOrder[i];
         vector<vector<NodeSharedPtr> > &tmp = elNodes[elType];
 
         for (j = 0; j < tmp.size(); ++j)
         {
             vector<int> tags;
             auto compIt = domainComposite.find(elType);
             if (compIt == domainComposite.end())
             {
                 tags.push_back(nComposite);
                 domainComposite[elType] = nComposite;
                 nComposite++;
             }
             else
             {
                 tags.push_back(compIt->second);
             }
 
             elMap[elIds[elType][j]] = reorderedId++;
 
             vector<NodeSharedPtr> nodeList = tmp[j];
 
             for (k = 0; k < nodeList.size(); ++k)
             {
                 auto testIns = m_mesh->m_vertexSet.insert(nodeList[k]);
 
                 if (!testIns.second)
                 {
                     nodeList[k] = *(testIns.first);
                 }
                 else
                 {
                     nodeList[k]->m_id = nodeCounter++;
                 }
             }
 
             // Create linear element
             ElmtConfig conf(elType, 1, false, false);
             ElementSharedPtr E = GetElementFactory().CreateInstance(
                 elType, conf, nodeList, tags);
             m_mesh->m_element[E->GetDim()].push_back(E);
         }
     }
 
     // -- Read in curved data.
     getline(m_mshFile, line);
     if (line.find("CURVE") == string::npos)
     {
         cerr << "Cannot find curved side data." << endl;
         abort();
     }
 
     // Read number of curves.
     getline(m_mshFile, line);
     s.clear();
     s.str(line);
     s >> nCurves;
 
     if (nCurves > 0)
     {
         string curveTag;
 
         for (i = 0; i < nCurves; ++i)
         {
             getline(m_mshFile, line);
             s.clear();
             s.str(line);
             s >> word;
 
             if (word == "File")
             {
                 // Next line contains filename and curve tag.
                 getline(m_mshFile, line);
                 s.clear();
                 s.str(line);
                 s >> word >> curveTag;
                 curveTags[curveTag] = make_pair(eFile, word);
             }
             else if (word == "Recon")
             {
                 // Next line contains curve tag.
                 getline(m_mshFile, line);
                 s.clear();
                 s.str(line);
                 s >> word >> curveTag;
                 curveTags[curveTag] = make_pair(eRecon, word);
             }
             else
             {
                 cerr << "Unsupported curve type " << word << endl;
                 abort();
             }
         }
 
         // Load high order surface information.
         LoadHOSurfaces();
 
         // Read in curve information. First line should contain number
         // of curved sides.
         getline(m_mshFile, line);
 
         if (line.find("side") == string::npos)
         {
             cerr << "Unable to read number of curved sides" << endl;
             abort();
         }
 
         int nCurvedSides;
         int faceId, elId;
 
         s.clear();
         s.str(line);
         s >> nCurvedSides;
 
         // Iterate over curved sides, and look up high-order surface
         // information in the HOSurfSet, then map this onto faces.
         for (i = 0; i < nCurvedSides; ++i)
         {
             getline(m_mshFile, line);
             s.clear();
             s.str(line);
             s >> faceId >> elId >> word;
             faceId--;
             elId                = elMap[elId - 1];
             ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][elId];
 
             int origFaceId = faceId;
 
             if (el->GetConf().m_e == LibUtilities::ePrism && faceId % 2 == 0)
             {
                 std::shared_ptr<Prism> p =
                     std::dynamic_pointer_cast<Prism>(el);
                 if (p->m_orientation == 1)
                 {
                     faceId = (faceId + 2) % 6;
                 }
                 else if (p->m_orientation == 2)
                 {
                     faceId = (faceId + 4) % 6;
                 }
             }
             else if (el->GetConf().m_e == LibUtilities::eTetrahedron)
             {
                 std::shared_ptr<Tetrahedron> t =
                     std::dynamic_pointer_cast<Tetrahedron>(el);
                 faceId = t->m_orientationMap[faceId];
             }
 
             auto it = curveTags.find(word);
             if (it == curveTags.end())
             {
                 cerr << "Unrecognised curve tag " << word << " in curved lines"
                      << endl;
                 abort();
             }
 
             if (it->second.first == eRecon)
             {
                 // Spherigon information: read in vertex normals.
                 vector<NodeSharedPtr> &tmp = el->GetFace(faceId)->m_vertexList;
                 vector<Node> n(tmp.size());
 
                 int offset = 0;
                 if (el->GetConf().m_e == LibUtilities::ePrism &&
                     faceId % 2 == 1)
                 {
                     offset =
                         std::dynamic_pointer_cast<Prism>(el)->m_orientation;
                 }
 
                 // Read x/y/z coordinates.
                 getline(m_mshFile, line);
                 s.clear();
                 s.str(line);
                 for (j = 0; j < tmp.size(); ++j)
                 {
                     s >> n[j].m_x;
                 }
 
                 getline(m_mshFile, line);
                 s.clear();
                 s.str(line);
                 for (j = 0; j < tmp.size(); ++j)
                 {
                     s >> n[j].m_y;
                 }
 
                 getline(m_mshFile, line);
                 s.clear();
                 s.str(line);
                 for (j = 0; j < tmp.size(); ++j)
                 {
                     s >> n[j].m_z;
                 }
 
                 for (j = 0; j < tmp.size(); ++j)
                 {
                     int id = tmp[(j + offset) % tmp.size()]->m_id;
                     auto vIt = m_mesh->m_vertexNormals.find(id);
 
                     if (vIt == m_mesh->m_vertexNormals.end())
                     {
                         m_mesh->m_vertexNormals[id] = n[j];
                     }
                 }
 
                 // Add edge/face to list of faces to apply spherigons
                 // to.
                 m_mesh->m_spherigonSurfs.insert(make_pair(elId, faceId));
             }
             else if (it->second.first == eFile)
             {
                 FaceSharedPtr f               = el->GetFace(faceId);
                 static int tetFaceVerts[4][3] = {
                     {0, 1, 2}, {0, 1, 3}, {1, 2, 3}, {0, 2, 3}};
 
                 vector<int> vertId(3);
                 s >> vertId[0] >> vertId[1] >> vertId[2];
 
                 // Prisms and tets may have been rotated by OrientPrism
                 // routine which reorders vertices. This block rotates
                 // vertex ids accordingly.
                 if (el->GetConf().m_e == LibUtilities::eTetrahedron)
                 {
                     std::shared_ptr<Tetrahedron> tet =
                         std::static_pointer_cast<Tetrahedron>(el);
                     vector<int> tmpVertId = vertId;
 
                     for (j = 0; j < 3; ++j)
                     {
                         int v =
                             tet->GetVertex(tet->m_origVertMap
                                                [tetFaceVerts[origFaceId][j]])
                                 ->m_id;
 
                         for (k = 0; k < 3; ++k)
                         {
                             int w = f->m_vertexList[k]->m_id;
                             if (v == w)
                             {
                                 vertId[k] = tmpVertId[j];
                                 break;
                             }
                         }
                     }
                 }
                 else if (el->GetConf().m_e == LibUtilities::ePrism)
                 {
                     std::shared_ptr<Prism> pr =
                         std::static_pointer_cast<Prism>(el);
                     if (pr->m_orientation == 1)
                     {
                         swap(vertId[2], vertId[1]);
                         swap(vertId[0], vertId[1]);
                     }
                     else if (pr->m_orientation == 2)
                     {
                         swap(vertId[0], vertId[1]);
                         swap(vertId[2], vertId[1]);
                     }
                 }
 
                 HOSurfSharedPtr hs =
                     std::shared_ptr<HOSurf>(new HOSurf(vertId));
                 // Find vertex combination in hoData.
                 auto hoIt = hoData[word].find(hs);
 
                 if (hoIt == hoData[word].end())
                 {
                     cerr << "Unable to find high-order surface data "
                          << "for element id " << elId + 1 << endl;
                     abort();
                 }
 
                 // Depending on order of vertices in rea file, surface
                 // information may need to be rotated or reflected.
                 HOSurfSharedPtr surf = *hoIt;
                 surf->Align(vertId);
 
                 // Finally, add high order data to appropriate face.
                 int Ntot = (*hoIt)->surfVerts.size();
                 int N    = ((int)sqrt(8.0 * Ntot + 1.0) - 1) / 2;
                 EdgeSharedPtr edge;
 
                 // Apply high-order map to convert face data to Nektar++
                 // ordering (vertices->edges->internal).
                 vector<NodeSharedPtr> tmpVerts = (*hoIt)->surfVerts;
                 for (j = 0; j < tmpVerts.size(); ++j)
                 {
                     (*hoIt)->surfVerts[hoMap[j]] = tmpVerts[j];
                 }
 
                 for (j = 0; j < tmpVerts.size(); ++j)
                 {
                     NodeSharedPtr a = (*hoIt)->surfVerts[j];
                 }
 
                 vector<int> faceVertIds(3);
                 faceVertIds[0] = f->m_vertexList[0]->m_id;
                 faceVertIds[1] = f->m_vertexList[1]->m_id;
                 faceVertIds[2] = f->m_vertexList[2]->m_id;
 
                 for (j = 0; j < f->m_edgeList.size(); ++j)
                 {
                     edge = f->m_edgeList[j];
 
                     // Skip over edges which have already been
                     // populated.
                     if (edge->m_edgeNodes.size() > 0)
                     {
                         continue;
                     }
 
                     // edge->m_edgeNodes.clear();
                     edge->m_curveType = LibUtilities::eGaussLobattoLegendre;
 
                     for (int k = 0; k < N - 2; ++k)
                     {
                         edge->m_edgeNodes.push_back(
                             (*hoIt)->surfVerts[3 + j * (N - 2) + k]);
                     }
 
                     // Nodal triangle data is always
                     // counter-clockwise. Add this check to reorder
                     // where necessary.
                     if (edge->m_n1->m_id != faceVertIds[j])
                     {
                         reverse(edge->m_edgeNodes.begin(),
                                 edge->m_edgeNodes.end());
                     }
                 }
 
                 // Insert interior face curvature.
                 f->m_curveType = LibUtilities::eNodalTriElec;
                 for (int j = 3 + 3 * (N - 2); j < Ntot; ++j)
                 {
                     f->m_faceNodes.push_back((*hoIt)->surfVerts[j]);
                 }
             }
         }
     }
 
     // -- Process fluid boundary conditions.
 
     // Define a fairly horrendous map: key is the condition ID, the
     // value is a vector of pairs of composites and element
     // types. Essentially this map takes conditions -> composites for
     // each element type.
     map<int, vector<pair<int, LibUtilities::ShapeType> > > surfaceCompMap;
 
     // Skip boundary conditions line.
     getline(m_mshFile, line);
     getline(m_mshFile, line);
 
     int nSurfaces = 0;
 
     while (true)
     {
         getline(m_mshFile, line);
 
         // Break out of loop at end of boundary conditions section.
         if (line.find("*") != string::npos || m_mshFile.eof() ||
             line.length() == 0)
         {
             break;
         }
 
         // Read boundary type, element ID and face ID.
         char bcType;
         int elId, faceId;
         s.clear();
         s.str(line);
         s >> bcType >> elId >> faceId;
         faceId--;
         elId = elMap[elId - 1];
 
         vector<string> vals;
         vector<ConditionType> type;
         ConditionSharedPtr c = MemoryManager<Condition>::AllocateSharedPtr();
 
         ElementSharedPtr elm = m_mesh->m_element[m_mesh->m_spaceDim][elId];
 
         // Ignore BCs for undefined edges/faces
         if ((elm->GetDim() == 2 && faceId >= elm->GetEdgeCount()) ||
             (elm->GetDim() == 3 && faceId >= elm->GetFaceCount()))
         {
             continue;
         }
 
         // First character on each line describes type of BC. Currently
         // only support V, W, and O. In this switch statement we
         // construct the quantities needed to search for the condition.
         switch (bcType)
         {
             // Wall boundary.
             case 'W':
             {
                 if (scalar)
                 {
                     vals.push_back("0");
                     type.push_back(eDirichlet);
                 }
                 else
                 {
                     for (i = 0; i < m_mesh->m_fields.size() - 1; ++i)
                     {
                         vals.push_back("0");
                         type.push_back(eDirichlet);
                     }
                     // Set high-order boundary condition for wall.
                     vals.push_back("0");
                     type.push_back(eHOPCondition);
                 }
                 break;
             }
 
             // Velocity boundary condition (either constant or dependent
             // upon x,y,z).
             case 'V':
             case 'v':
             {
                 if (scalar)
                 {
                     getline(m_mshFile, line);
                     size_t p = line.find_first_of('=');
                     vals.push_back(
                         boost::algorithm::trim_copy(line.substr(p + 1)));
                     type.push_back(eDirichlet);
                 }
                 else
                 {
                     for (i = 0; i < m_mesh->m_fields.size() - 1; ++i)
                     {
                         getline(m_mshFile, line);
                         size_t p = line.find_first_of('=');
                         vals.push_back(
                             boost::algorithm::trim_copy(line.substr(p + 1)));
                         type.push_back(eDirichlet);
                     }
                     // Set high-order boundary condition for Dirichlet
                     // condition.
                     vals.push_back("0");
                     type.push_back(eHOPCondition);
                 }
                 break;
             }
 
             // Natural outflow condition (default value = 0.0?)
             case 'O':
             {
                 if (scalar)
                 {
                     vals.push_back("0");
                     type.push_back(eNeumann);
                 }
                 else
                 {
                     for (i = 0; i < m_mesh->m_fields.size(); ++i)
                     {
                         vals.push_back("0");
                         type.push_back(eNeumann);
                     }
                     // Set zero Dirichlet condition for outflow.
                     type[m_mesh->m_fields.size() - 1] = eDirichlet;
                 }
                 break;
             }
 
             // Ignore unsupported BCs
             case 'P':
             case 'E':
             case 'F':
             case 'f':
             case 'N':
                 continue;
                 break;
 
             default:
                 cerr << "Unknown boundary condition type " << line[0] << endl;
                 abort();
         }
 
         // Populate condition information.
         c->field = m_mesh->m_fields;
         c->type  = type;
         c->value = vals;
 
         // Now attempt to find this boundary condition inside
         // m_mesh->condition. This is currently a linear search and should
         // probably be made faster!
         bool found = false;
         auto it = m_mesh->m_condition.begin();
         for (; it != m_mesh->m_condition.end(); ++it)
         {
             if (c == it->second)
             {
                 found = true;
                 break;
             }
         }
 
         int compTag, conditionId;
         ElementSharedPtr surfEl;
 
         // Create element for face (3D) or segment (2D). At the moment
         // this is a bit of a hack since high-order nodes are not
         // copied, so some output modules (e.g. Gmsh) will not output
         // correctly.
         if (elm->GetDim() == 3)
         {
             // 3D elements may have been reoriented, so face IDs will
             // change.
             if (elm->GetConf().m_e == LibUtilities::ePrism && faceId % 2 == 0)
             {
                 std::shared_ptr<Prism> p =
                     std::dynamic_pointer_cast<Prism>(elm);
                 if (p->m_orientation == 1)
                 {
                     faceId = (faceId + 2) % 6;
                 }
                 else if (p->m_orientation == 2)
                 {
                     faceId = (faceId + 4) % 6;
                 }
             }
             else if (elm->GetConf().m_e == LibUtilities::eTetrahedron)
             {
                 std::shared_ptr<Tetrahedron> t =
                     std::dynamic_pointer_cast<Tetrahedron>(elm);
                 faceId = t->m_orientationMap[faceId];
             }
 
             FaceSharedPtr f = elm->GetFace(faceId);
             bool tri        = f->m_vertexList.size() == 3;
 
             vector<NodeSharedPtr> nodeList;
             nodeList.insert(nodeList.begin(),
                             f->m_vertexList.begin(),
                             f->m_vertexList.end());
 
             vector<int> tags;
 
             LibUtilities::ShapeType seg =
                 tri ? LibUtilities::eTriangle : LibUtilities::eQuadrilateral;
             ElmtConfig conf(
                 seg, 1, true, true, false, LibUtilities::eGaussLobattoLegendre);
             surfEl =
                 GetElementFactory().CreateInstance(seg, conf, nodeList, tags);
 
             // Copy high-order surface information from edges.
             for (int i = 0; i < f->m_vertexList.size(); ++i)
             {
                 surfEl->GetEdge(i)->m_edgeNodes = f->m_edgeList[i]->m_edgeNodes;
                 surfEl->GetEdge(i)->m_curveType = f->m_edgeList[i]->m_curveType;
             }
         }
         else if (faceId < elm->GetEdgeCount())
         {
             EdgeSharedPtr f = elm->GetEdge(faceId);
 
             vector<NodeSharedPtr> nodeList;
             nodeList.push_back(f->m_n1);
             nodeList.push_back(f->m_n2);
 
             vector<int> tags;
 
             ElmtConfig conf(LibUtilities::eSegment,
                             1,
                             true,
                             true,
                             false,
                             LibUtilities::eGaussLobattoLegendre);
             surfEl = GetElementFactory().CreateInstance(
                 LibUtilities::eSegment, conf, nodeList, tags);
         }
 
         if (!surfEl)
         {
             continue;
         }
 
         LibUtilities::ShapeType surfElType = surfEl->GetConf().m_e;
 
         if (!found)
         {
             // If condition does not already exist, add to condition
             // list, create new composite tag and put inside
             // surfaceCompMap.
             conditionId = m_mesh->m_condition.size();
             compTag = nComposite;
             c->m_composite.push_back(compTag);
             m_mesh->m_condition[conditionId] = c;
 
             surfaceCompMap[conditionId].push_back(
                 pair<int, LibUtilities::ShapeType>(nComposite, surfElType));
 
             nComposite++;
         }
         else
         {
             // Otherwise find existing composite inside surfaceCompMap.
             auto it2 = surfaceCompMap.find(it->first);
 
             found = false;
             if (it2 == surfaceCompMap.end())
             {
                 // This should never happen!
                 cerr << "Unable to find condition!" << endl;
                 abort();
             }
 
             for (j = 0; j < it2->second.size(); ++j)
             {
                 pair<int, LibUtilities::ShapeType> tmp = it2->second[j];
                 if (tmp.second == surfElType)
                 {
                     found   = true;
                     compTag = tmp.first;
                     break;
                 }
             }
 
             // If no pairs were found, then this condition contains
             // multiple element types (i.e. both triangles and
             // quads). Create another composite for the new shape type
             // and insert into the map.
             if (!found)
             {
                 it2->second.push_back(
                     pair<int, LibUtilities::ShapeType>(nComposite, surfElType));
                 compTag = nComposite;
                 m_mesh->m_condition[it->first]->m_composite.push_back(compTag);
                 nComposite++;
             }
 
             conditionId = it->first;
         }
 
         // Insert composite tag into element and insert element into
         // mesh.
         vector<int> existingTags = surfEl->GetTagList();
 
         existingTags.insert(existingTags.begin(), compTag);
         surfEl->SetTagList(existingTags);
         surfEl->SetId(nSurfaces);
 
         m_mesh->m_element[surfEl->GetDim()].push_back(surfEl);
         nSurfaces++;
     }
 
     m_mshFile.reset();
 
     // -- Process rest of mesh.
     ProcessEdges();
     ProcessFaces();
     ProcessElements();
     ProcessComposites();
 }

Member Data Documentation

◆ className

ModuleKey Nektar::Utilities::InputNek::className

static

Initial value:

= GetModuleFactory().RegisterCreatorFunction(

ModuleKey(eInputModule, "rea"), InputNek::create, "Reads Nektar rea file.")

ModuleKey for class.

Definition at line 72 of file InputNek.h.

◆ curveTags

std::map<std::string, std::pair<NekCurve, std::string> > Nektar::Utilities::InputNek::curveTags

private

Maps a curve tag to a filename containing surface information.

Definition at line 81 of file InputNek.h.

Referenced by LoadHOSurfaces(), and Process().

◆ hoData

std::map<std::string, NekMeshUtils::HOSurfSet> Nektar::Utilities::InputNek::hoData

private

Maps a curve tag to the high-order surface data for that tag.

Definition at line 86 of file InputNek.h.

Referenced by LoadHOSurfaces(), and Process().

◆ hoMap

std::map<int, int> Nektar::Utilities::InputNek::hoMap

private

Maps ordering of hsf standard element to Nektar++ ordering.

Definition at line 91 of file InputNek.h.

Referenced by LoadHOSurfaces(), and Process().

Public Member Functions

Static Public Member Functions

Static Public Attributes

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ InputNek()

◆ ~InputNek()

Member Function Documentation

◆ create()

◆ GetNnodes()

◆ LoadHOSurfaces()

◆ Process()

Member Data Documentation

◆ className

◆ curveTags

◆ hoData

◆ hoMap