Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Static Public Member Functions | Static Public Attributes | Private Member Functions | Private Attributes | List of all members
Nektar::Utilities::InputNek Class Reference

#include <InputNek.h>

Inheritance diagram for Nektar::Utilities::InputNek:
Inheritance graph
[legend]
Collaboration diagram for Nektar::Utilities::InputNek:
Collaboration graph
[legend]

Public Member Functions

 InputNek (MeshSharedPtr p_m)
virtual ~InputNek ()
virtual void Process ()
 Processes Nektar file format.
- Public Member Functions inherited from Nektar::Utilities::InputModule
 InputModule (FieldSharedPtr p_m)
void AddFile (string fileType, string fileName)
 InputModule (MeshSharedPtr p_m)
void OpenStream ()
 Open a file for input.
- Public Member Functions inherited from Nektar::Utilities::Module
 Module (FieldSharedPtr p_f)
virtual void Process (po::variables_map &vm)=0
void RegisterConfig (string key, string value)
 Register a configuration option with a module.
void PrintConfig ()
 Print out all configuration options for a module.
void SetDefaults ()
 Sets default configuration options for those which have not been set.
bool GetRequireEquiSpaced (void)
void SetRequireEquiSpaced (bool pVal)
void EvaluateTriFieldAtEquiSpacedPts (LocalRegions::ExpansionSharedPtr &exp, const Array< OneD, const NekDouble > &infield, Array< OneD, NekDouble > &outfield)
 Module (MeshSharedPtr p_m)
void RegisterConfig (string key, string value)
void PrintConfig ()
void SetDefaults ()
MeshSharedPtr GetMesh ()
virtual void ProcessVertices ()
 Extract element vertices.

Static Public Member Functions

static ModuleSharedPtr create (MeshSharedPtr m)
 Creates an instance of this class.

Static Public Attributes

static ModuleKey className
 ModuleKey for class.

Private Member Functions

void LoadHOSurfaces ()
int GetNnodes (LibUtilities::ShapeType elType)

Private Attributes

std::map< string, pair
< NekCurve, string > > 		curveTags
std::map< string, HOSurfSet > hoData
std::map< int, int > hoMap

Additional Inherited Members

- Protected Member Functions inherited from Nektar::Utilities::InputModule
void PrintSummary ()
 Print summary of elements.
void PrintSummary ()
 Print summary of elements.
- Protected Attributes inherited from Nektar::Utilities::InputModule
set< string > m_allowedFiles
std::ifstream m_mshFile
 Input stream.

Detailed Description

Converter class for Nektar session files.

Definition at line 83 of file InputNek.h.

Constructor & Destructor Documentation

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

Definition at line 57 of file InputNek.cpp.

: InputModule(m)
{
}
Nektar::Utilities::InputNek::~InputNek ( )
virtual

Definition at line 62 of file InputNek.cpp.

{
}

Member Function Documentation

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

Creates an instance of this class.

Definition at line 91 of file InputNek.h.

{
return MemoryManager<InputNek>::AllocateSharedPtr(m);
}
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 1057 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;
}
void Nektar::Utilities::InputNek::LoadHOSurfaces ( )
private

Load high order surface information from hsf file.

Definition at line 895 of file InputNek.cpp.

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

Referenced by Process().

{
map<string, pair<NekCurve, string> >::iterator it;
int nodeId = m_mesh->GetNumEntities();
for (it = curveTags.begin(); it != curveTags.end(); ++it)
{
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();
}
}
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
pFilenameFilename of Nektar session file to read.

Implements Nektar::Utilities::Module.

Definition at line 77 of file InputNek.cpp.

References curveTags, Nektar::Utilities::eDirichlet, Nektar::Utilities::eFile, Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eHexahedron, Nektar::Utilities::eHOPCondition, Nektar::Utilities::eNeumann, Nektar::LibUtilities::eNodalTriElec, Nektar::LibUtilities::ePrism, Nektar::LibUtilities::ePyramid, Nektar::LibUtilities::eQuadrilateral, Nektar::Utilities::eRecon, Nektar::LibUtilities::eSegment, Nektar::LibUtilities::eTetrahedron, Nektar::LibUtilities::eTriangle, Nektar::Utilities::GetElementFactory(), GetNnodes(), Nektar::Utilities::Element::GetVertex(), hoData, hoMap, Nektar::iterator, LoadHOSurfaces(), Nektar::Utilities::Module::m_mesh, Nektar::Utilities::InputModule::m_mshFile, Nektar::Utilities::Tetrahedron::m_orientationMap, Nektar::Utilities::InputModule::OpenStream(), Nektar::Utilities::Module::ProcessComposites(), Nektar::Utilities::Module::ProcessEdges(), Nektar::Utilities::Module::ProcessElements(), and Nektar::Utilities::Module::ProcessFaces().

{
// 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;
boost::unordered_map<int,int> elMap;
vector<LibUtilities::ShapeType> elmOrder;
// 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");
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 = boost::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;
map<LibUtilities::ShapeType,int>::iterator 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)
{
pair<NodeSet::iterator, bool> 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;
map<string,pair<NekCurve, string> >::iterator it;
HOSurfSet::iterator hoIt;
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)
{
boost::shared_ptr<Prism> p =
boost::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)
{
boost::shared_ptr<Tetrahedron> t =
boost::dynamic_pointer_cast<Tetrahedron>(el);
faceId = t->m_orientationMap[faceId];
}
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 = boost::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;
boost::unordered_map<int, Node>::iterator 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)
{
boost::shared_ptr<Tetrahedron> tet =
boost::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)
{
boost::shared_ptr<Prism> pr =
boost::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]);
}
}
// Find vertex combination in hoData.
hoIt = hoData[word].find(HOSurfSharedPtr(
new HOSurf(vertId)));
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();
// 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':
{
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':
{
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':
{
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!
ConditionMap::iterator it;
bool found = false;
for (it = m_mesh->m_condition.begin(); it != m_mesh->m_condition.end(); ++it)
{
if (c == it->second)
{
found = true;
break;
}
}
int compTag, conditionId;
ElementSharedPtr elm = m_mesh->m_element[m_mesh->m_spaceDim][elId];
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)
{
boost::shared_ptr<Prism> p =
boost::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)
{
boost::shared_ptr<Tetrahedron> t =
boost::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
{
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 (!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,surfEl->GetConf().m_e));
nComposite++;
}
else
{
// Otherwise find existing composite inside surfaceCompMap.
map<int,vector<pair<int,LibUtilities::ShapeType> > >::iterator it2;
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 == surfEl->GetConf().m_e)
{
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,surfEl->GetConf().m_e));
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.close();
// -- Process rest of mesh.
ProcessEdges ();
ProcessFaces ();
ProcessElements ();
ProcessComposites();
}

Member Data Documentation

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 95 of file InputNek.h.

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

Maps a curve tag to a filename containing surface information.

Definition at line 104 of file InputNek.h.

Referenced by LoadHOSurfaces(), and Process().

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

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

Definition at line 109 of file InputNek.h.

Referenced by LoadHOSurfaces(), and Process().

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

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

Definition at line 114 of file InputNek.h.

Referenced by LoadHOSurfaces(), and Process().

Generated on Wed Mar 4 2015 16:33:50 for Nektar++ by   doxygen 1.8.1.2