Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 (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)
 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 60 of file InputNek.h.

Constructor & Destructor Documentation

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

Definition at line 60 of file InputNek.cpp.

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

60  : InputModule(m)
61 {
62  m_config["scalar"] = ConfigOption(
63  true, "0", "If defined then assume input rea is for scalar problem");
64 }
Represents a command-line configuration option.
std::map< std::string, ConfigOption > m_config
List of configuration values.
NEKMESHUTILS_EXPORT InputModule(MeshSharedPtr p_m)
Nektar::Utilities::InputNek::~InputNek ( )
virtual

Definition at line 66 of file InputNek.cpp.

67 {
68 }

Member Function Documentation

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

Creates an instance of this class.

Definition at line 68 of file InputNek.h.

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

69  {
71  }
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
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 1124 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().

1125 {
1126  int nNodes = 0;
1127 
1128  switch (InputNekEntity)
1129  {
1130  case LibUtilities::ePoint:
1131  nNodes = 1;
1132  break;
1134  nNodes = 2;
1135  break;
1137  nNodes = 3;
1138  break;
1140  nNodes = 4;
1141  break;
1143  nNodes = 4;
1144  break;
1146  nNodes = 5;
1147  break;
1148  case LibUtilities::ePrism:
1149  nNodes = 6;
1150  break;
1152  nNodes = 8;
1153  break;
1154  default:
1155  cerr << "unknown Nektar element type" << endl;
1156  }
1157 
1158  return nNodes;
1159 }
void Nektar::Utilities::InputNek::LoadHOSurfaces ( )
private

Load high order surface information from hsf file.

Definition at line 959 of file InputNek.cpp.

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

Referenced by Process().

960 {
961  map<string, pair<NekCurve, string> >::iterator it;
962  int nodeId = m_mesh->GetNumEntities();
963 
964  for (it = curveTags.begin(); it != curveTags.end(); ++it)
965  {
966  ifstream hsf;
967  string line, fileName = it->second.second;
968  size_t pos;
969  int N, Nface, dot;
970 
971  if (it->second.first != eFile)
972  {
973  continue;
974  }
975 
976  // Replace fro extension with hsf.
977  dot = fileName.find_last_of('.');
978  fileName = fileName.substr(0, dot);
979  fileName += ".hsf";
980 
981  // Open hsf file.
982  hsf.open(fileName.c_str());
983  if (!hsf.is_open())
984  {
985  cerr << "Could not open surface file " << fileName << endl;
986  abort();
987  }
988 
989  // Read in header line; determine element order, number of faces
990  // from this line.
991  getline(hsf, line);
992  pos = line.find("=");
993  if (pos == string::npos)
994  {
995  cerr << "hsf header error: cannot read number of "
996  << "nodal points." << endl;
997  abort();
998  }
999  line = line.substr(pos + 1);
1000  stringstream ss(line);
1001  ss >> N;
1002 
1003  pos = line.find("=");
1004  if (pos == string::npos)
1005  {
1006  cerr << "hsf header error: cannot read number of "
1007  << "faces." << endl;
1008  abort();
1009  }
1010  line = line.substr(pos + 1);
1011  ss.clear();
1012  ss.str(line);
1013  ss >> Nface;
1014 
1015  int Ntot = N * (N + 1) / 2;
1016 
1017  // Skip a line, then read in r,s positions inside the next
1018  // comments.
1019  Array<OneD, NekDouble> r(Ntot), s(Ntot);
1020  getline(hsf, line);
1021 
1022  for (int i = 0; i < 2; ++i)
1023  {
1024  string word;
1025 
1026  getline(hsf, line);
1027  ss.clear();
1028  ss.str(line);
1029  ss >> word;
1030 
1031  if (word != "#")
1032  {
1033  cerr << "hsf header error: cannot read in "
1034  << "r/s points" << endl;
1035  abort();
1036  }
1037 
1038  for (int j = 0; j < Ntot; ++j)
1039  {
1040  ss >> (i == 0 ? r[j] : s[j]);
1041  }
1042  }
1043 
1044  // Generate electrostatic points so that re-mapping array can
1045  // be constructed.
1046  Array<OneD, NekDouble> rp(Ntot), sp(Ntot);
1047  LibUtilities::PointsKey elec(N, LibUtilities::eNodalTriElec);
1048  LibUtilities::PointsManager()[elec]->GetPoints(rp, sp);
1049 
1050  // Expensively construct remapping array nodemap. This will
1051  // map nodal ordering from hsf order to Nektar++ ordering
1052  // (i.e. vertices followed by edges followed by interior
1053  // points.)
1054  for (int i = 0; i < Ntot; ++i)
1055  {
1056  for (int j = 0; j < Ntot; ++j)
1057  {
1058  if (fabs(r[i] - rp[j]) < 1e-5 && fabs(s[i] - sp[j]) < 1e-5)
1059  {
1060  hoMap[i] = j;
1061  break;
1062  }
1063  }
1064  }
1065 
1066  // Skip variables line
1067  getline(hsf, line);
1068 
1069  // Read in nodal points for each face.
1070  map<int, vector<NodeSharedPtr> > faceMap;
1071  for (int i = 0; i < Nface; ++i)
1072  {
1073  getline(hsf, line);
1074  vector<NodeSharedPtr> faceNodes(Ntot);
1075  for (int j = 0; j < Ntot; ++j, ++nodeId)
1076  {
1077  double x, y, z;
1078  getline(hsf, line);
1079  ss.clear();
1080  ss.str(line);
1081  ss >> x >> y >> z;
1082  faceNodes[j] = NodeSharedPtr(new Node(nodeId, x, y, z));
1083  }
1084  // Skip over tecplot connectivity information.
1085  for (int j = 0; j < (N - 1) * (N - 1); ++j)
1086  {
1087  getline(hsf, line);
1088  }
1089  faceMap[i] = faceNodes;
1090  }
1091 
1092  // Finally, read in connectivity information to set up after
1093  // reading rea file.
1094  getline(hsf, line);
1095  for (int i = 0; i < Nface; ++i)
1096  {
1097  string tmp;
1098  int fid;
1099  vector<int> nodeIds(3);
1100 
1101  getline(hsf, line);
1102  ss.clear();
1103  ss.str(line);
1104  ss >> tmp >> fid >> nodeIds[0] >> nodeIds[1] >> nodeIds[2];
1105 
1106  if (tmp != "#")
1107  {
1108  cerr << "Unable to read hsf connectivity information." << endl;
1109  abort();
1110  }
1111 
1112  hoData[it->first].insert(
1113  HOSurfSharedPtr(new HOSurf(nodeIds, faceMap[i])));
1114  }
1115 
1116  hsf.close();
1117  }
1118 }
boost::shared_ptr< HOSurf > HOSurfSharedPtr
Definition: HOAlignment.h:181
std::map< int, int > hoMap
Definition: InputNek.h:92
std::map< std::string, NekMeshUtils::HOSurfSet > hoData
Definition: InputNek.h:87
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
PointsManagerT & PointsManager(void)
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
std::map< std::string, std::pair< NekCurve, std::string > > curveTags
Definition: InputNek.h:82
HOTriangle< NodeSharedPtr > HOSurf
Definition: HOAlignment.h:180
2D Nodal Electrostatic Points on a Triangle
Definition: PointsType.h:70
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::NekMeshUtils::Module.

Definition at line 79 of file InputNek.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), curveTags, Nektar::SpatialDomains::eDirichlet, Nektar::Utilities::eFile, Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eHexahedron, Nektar::NekMeshUtils::eHOPCondition, Nektar::SpatialDomains::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::NekMeshUtils::GetElementFactory(), GetNnodes(), Nektar::NekMeshUtils::Element::GetVertex(), hoData, hoMap, Nektar::iterator, LoadHOSurfaces(), Nektar::NekMeshUtils::Module::m_config, Nektar::NekMeshUtils::Module::m_mesh, Nektar::NekMeshUtils::InputModule::m_mshFile, Nektar::NekMeshUtils::Tetrahedron::m_orientationMap, class_topology::Node, Nektar::NekMeshUtils::InputModule::OpenStream(), CellMLToNektar.cellml_metadata::p, Nektar::NekMeshUtils::Module::ProcessComposites(), Nektar::NekMeshUtils::Module::ProcessEdges(), Nektar::NekMeshUtils::Module::ProcessElements(), and Nektar::NekMeshUtils::Module::ProcessFaces().

80 {
81  // Open the file stream.
82  OpenStream();
83 
84  string line, word;
85  int nParam, nElements, nCurves;
86  int i, j, k, nodeCounter = 0;
87  int nComposite = 0;
89  double vertex[3][8];
90  map<LibUtilities::ShapeType, int> domainComposite;
91  map<LibUtilities::ShapeType, vector<vector<NodeSharedPtr> > > elNodes;
92  map<LibUtilities::ShapeType, vector<int> > elIds;
93  boost::unordered_map<int, int> elMap;
94  vector<LibUtilities::ShapeType> elmOrder;
95 
96  bool scalar = m_config["scalar"].as<bool>();
97 
98  // Set up vector of processing orders.
99  elmOrder.push_back(LibUtilities::eSegment);
100  elmOrder.push_back(LibUtilities::eTriangle);
101  elmOrder.push_back(LibUtilities::eQuadrilateral);
102  elmOrder.push_back(LibUtilities::ePrism);
103  elmOrder.push_back(LibUtilities::ePyramid);
104  elmOrder.push_back(LibUtilities::eTetrahedron);
105  elmOrder.push_back(LibUtilities::eHexahedron);
106 
107  m_mesh->m_expDim = 0;
108  m_mesh->m_spaceDim = 0;
109 
110  if (m_mesh->m_verbose)
111  {
112  cout << "InputNek: Start reading file..." << endl;
113  }
114 
115  // -- Read in parameters.
116 
117  // Ignore first 3 lines. 4th line contains number of parameters.
118  for (i = 0; i < 4; ++i)
119  {
120  getline(m_mshFile, line);
121  }
122 
123  stringstream s(line);
124  s >> nParam;
125 
126  for (i = 0; i < nParam; ++i)
127  {
128  string tmp1, tmp2;
129  getline(m_mshFile, line);
130  s.str(line);
131  s >> tmp1 >> tmp2;
132  }
133 
134  // -- Read in passive scalars (ignore)
135  getline(m_mshFile, line);
136  s.clear();
137  s.str(line);
138  s >> j;
139  for (i = 0; i < j; ++i)
140  {
141  getline(m_mshFile, line);
142  }
143 
144  // -- Read in logical switches (ignore)
145  getline(m_mshFile, line);
146  s.clear();
147  s.str(line);
148  s >> j;
149  for (i = 0; i < j; ++i)
150  {
151  getline(m_mshFile, line);
152  }
153 
154  // -- Read in mesh data.
155 
156  // First hunt for MESH tag
157  bool foundMesh = false;
158  while (!m_mshFile.eof())
159  {
160  getline(m_mshFile, line);
161  if (line.find("MESH") != string::npos)
162  {
163  foundMesh = true;
164  break;
165  }
166  }
167 
168  if (!foundMesh)
169  {
170  cerr << "Couldn't find MESH tag inside file." << endl;
171  abort();
172  }
173 
174  // Now read in number of elements and space dimension.
175  getline(m_mshFile, line);
176  s.clear();
177  s.str(line);
178  s >> nElements >> m_mesh->m_expDim;
179  m_mesh->m_spaceDim = m_mesh->m_expDim;
180 
181  // Set up field names.
182  m_mesh->m_fields.push_back("u");
183 
184  if (!scalar)
185  {
186  m_mesh->m_fields.push_back("v");
187  if (m_mesh->m_spaceDim > 2)
188  {
189  m_mesh->m_fields.push_back("w");
190  }
191  m_mesh->m_fields.push_back("p");
192  }
193 
194  // Loop over and create elements.
195  for (i = 0; i < nElements; ++i)
196  {
197  getline(m_mshFile, line);
198 
199  if (m_mesh->m_expDim == 2)
200  {
201  if (line.find("Qua") != string::npos ||
202  line.find("qua") != string::npos)
203  {
205  }
206  else
207  {
208  // Default element type in 2D is triangle.
209  elType = LibUtilities::eTriangle;
210  }
211  }
212  else
213  {
214  if (line.find("Tet") != string::npos ||
215  line.find("tet") != string::npos)
216  {
218  }
219  else if (line.find("Hex") != string::npos ||
220  line.find("hex") != string::npos)
221  {
222  elType = LibUtilities::eHexahedron;
223  }
224  else if (line.find("Prism") != string::npos ||
225  line.find("prism") != string::npos)
226  {
227  elType = LibUtilities::ePrism;
228  }
229  else if (line.find("Pyr") != string::npos ||
230  line.find("pyr") != string::npos)
231  {
232  elType = LibUtilities::ePyramid;
233  }
234  else if (line.find("Qua") != string::npos ||
235  line.find("qua") != string::npos)
236  {
238  }
239  else
240  {
241  // Default element type in 2D is tetrahedron.
243  }
244  }
245 
246  // Read in number of vertices for element type.
247  const int nNodes = GetNnodes(elType);
248 
249  for (j = 0; j < m_mesh->m_expDim; ++j)
250  {
251  getline(m_mshFile, line);
252  s.clear();
253  s.str(line);
254  for (k = 0; k < nNodes; ++k)
255  {
256  s >> vertex[j][k];
257  }
258  }
259 
260  // Zero co-ordinates bigger than expansion dimension.
261  for (j = m_mesh->m_expDim; j < 3; ++j)
262  {
263  for (k = 0; k < nNodes; ++k)
264  {
265  vertex[j][k] = 0.0;
266  }
267  }
268 
269  // Nektar meshes do not contain a unique list of nodes, so this
270  // block constructs a unique set so that elements can be created
271  // with unique nodes.
272  vector<NodeSharedPtr> nodeList;
273  for (k = 0; k < nNodes; ++k)
274  {
275  NodeSharedPtr n = boost::shared_ptr<Node>(new Node(
276  nodeCounter++, vertex[0][k], vertex[1][k], vertex[2][k]));
277  nodeList.push_back(n);
278  }
279 
280  elNodes[elType].push_back(nodeList);
281  elIds[elType].push_back(i);
282  }
283 
284  int reorderedId = 0;
285  nodeCounter = 0;
286 
287  for (i = 0; i < elmOrder.size(); ++i)
288  {
289  LibUtilities::ShapeType elType = elmOrder[i];
290  vector<vector<NodeSharedPtr> > &tmp = elNodes[elType];
291 
292  for (j = 0; j < tmp.size(); ++j)
293  {
294  vector<int> tags;
296  domainComposite.find(elType);
297  if (compIt == domainComposite.end())
298  {
299  tags.push_back(nComposite);
300  domainComposite[elType] = nComposite;
301  nComposite++;
302  }
303  else
304  {
305  tags.push_back(compIt->second);
306  }
307 
308  elMap[elIds[elType][j]] = reorderedId++;
309 
310  vector<NodeSharedPtr> nodeList = tmp[j];
311 
312  for (k = 0; k < nodeList.size(); ++k)
313  {
314  pair<NodeSet::iterator, bool> testIns =
315  m_mesh->m_vertexSet.insert(nodeList[k]);
316 
317  if (!testIns.second)
318  {
319  nodeList[k] = *(testIns.first);
320  }
321  else
322  {
323  nodeList[k]->m_id = nodeCounter++;
324  }
325  }
326 
327  // Create linear element
328  ElmtConfig conf(elType, 1, false, false);
330  elType, conf, nodeList, tags);
331  m_mesh->m_element[E->GetDim()].push_back(E);
332  }
333  }
334 
335  // -- Read in curved data.
336  getline(m_mshFile, line);
337  if (line.find("CURVE") == string::npos)
338  {
339  cerr << "Cannot find curved side data." << endl;
340  abort();
341  }
342 
343  // Read number of curves.
344  getline(m_mshFile, line);
345  s.clear();
346  s.str(line);
347  s >> nCurves;
348 
349  if (nCurves > 0)
350  {
351  string curveTag;
352 
353  for (i = 0; i < nCurves; ++i)
354  {
355  getline(m_mshFile, line);
356  s.clear();
357  s.str(line);
358  s >> word;
359 
360  if (word == "File")
361  {
362  // Next line contains filename and curve tag.
363  getline(m_mshFile, line);
364  s.clear();
365  s.str(line);
366  s >> word >> curveTag;
367  curveTags[curveTag] = make_pair(eFile, word);
368  }
369  else if (word == "Recon")
370  {
371  // Next line contains curve tag.
372  getline(m_mshFile, line);
373  s.clear();
374  s.str(line);
375  s >> word >> curveTag;
376  curveTags[curveTag] = make_pair(eRecon, word);
377  }
378  else
379  {
380  cerr << "Unsupported curve type " << word << endl;
381  abort();
382  }
383  }
384 
385  // Load high order surface information.
386  LoadHOSurfaces();
387 
388  // Read in curve information. First line should contain number
389  // of curved sides.
390  getline(m_mshFile, line);
391 
392  if (line.find("side") == string::npos)
393  {
394  cerr << "Unable to read number of curved sides" << endl;
395  abort();
396  }
397 
398  int nCurvedSides;
399  int faceId, elId;
400  map<string, pair<NekCurve, string> >::iterator it;
401  HOSurfSet::iterator hoIt;
402 
403  s.clear();
404  s.str(line);
405  s >> nCurvedSides;
406 
407  // Iterate over curved sides, and look up high-order surface
408  // information in the HOSurfSet, then map this onto faces.
409  for (i = 0; i < nCurvedSides; ++i)
410  {
411  getline(m_mshFile, line);
412  s.clear();
413  s.str(line);
414  s >> faceId >> elId >> word;
415  faceId--;
416  elId = elMap[elId - 1];
417  ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][elId];
418 
419  int origFaceId = faceId;
420 
421  if (el->GetConf().m_e == LibUtilities::ePrism && faceId % 2 == 0)
422  {
423  boost::shared_ptr<Prism> p =
424  boost::dynamic_pointer_cast<Prism>(el);
425  if (p->m_orientation == 1)
426  {
427  faceId = (faceId + 2) % 6;
428  }
429  else if (p->m_orientation == 2)
430  {
431  faceId = (faceId + 4) % 6;
432  }
433  }
434  else if (el->GetConf().m_e == LibUtilities::eTetrahedron)
435  {
436  boost::shared_ptr<Tetrahedron> t =
437  boost::dynamic_pointer_cast<Tetrahedron>(el);
438  faceId = t->m_orientationMap[faceId];
439  }
440 
441  it = curveTags.find(word);
442  if (it == curveTags.end())
443  {
444  cerr << "Unrecognised curve tag " << word << " in curved lines"
445  << endl;
446  abort();
447  }
448 
449  if (it->second.first == eRecon)
450  {
451  // Spherigon information: read in vertex normals.
452  vector<NodeSharedPtr> &tmp = el->GetFace(faceId)->m_vertexList;
453  vector<Node> n(tmp.size());
454 
455  int offset = 0;
456  if (el->GetConf().m_e == LibUtilities::ePrism &&
457  faceId % 2 == 1)
458  {
459  offset =
460  boost::dynamic_pointer_cast<Prism>(el)->m_orientation;
461  }
462 
463  // Read x/y/z coordinates.
464  getline(m_mshFile, line);
465  s.clear();
466  s.str(line);
467  for (j = 0; j < tmp.size(); ++j)
468  {
469  s >> n[j].m_x;
470  }
471 
472  getline(m_mshFile, line);
473  s.clear();
474  s.str(line);
475  for (j = 0; j < tmp.size(); ++j)
476  {
477  s >> n[j].m_y;
478  }
479 
480  getline(m_mshFile, line);
481  s.clear();
482  s.str(line);
483  for (j = 0; j < tmp.size(); ++j)
484  {
485  s >> n[j].m_z;
486  }
487 
488  for (j = 0; j < tmp.size(); ++j)
489  {
490  int id = tmp[(j + offset) % tmp.size()]->m_id;
492  m_mesh->m_vertexNormals.find(id);
493 
494  if (vIt == m_mesh->m_vertexNormals.end())
495  {
496  m_mesh->m_vertexNormals[id] = n[j];
497  }
498  }
499 
500  // Add edge/face to list of faces to apply spherigons
501  // to.
502  m_mesh->m_spherigonSurfs.insert(make_pair(elId, faceId));
503  }
504  else if (it->second.first == eFile)
505  {
506  FaceSharedPtr f = el->GetFace(faceId);
507  static int tetFaceVerts[4][3] = {
508  {0, 1, 2}, {0, 1, 3}, {1, 2, 3}, {0, 2, 3}};
509 
510  vector<int> vertId(3);
511  s >> vertId[0] >> vertId[1] >> vertId[2];
512 
513  // Prisms and tets may have been rotated by OrientPrism
514  // routine which reorders vertices. This block rotates
515  // vertex ids accordingly.
516  if (el->GetConf().m_e == LibUtilities::eTetrahedron)
517  {
518  boost::shared_ptr<Tetrahedron> tet =
519  boost::static_pointer_cast<Tetrahedron>(el);
520  vector<int> tmpVertId = vertId;
521 
522  for (j = 0; j < 3; ++j)
523  {
524  int v =
525  tet->GetVertex(tet->m_origVertMap
526  [tetFaceVerts[origFaceId][j]])
527  ->m_id;
528 
529  for (k = 0; k < 3; ++k)
530  {
531  int w = f->m_vertexList[k]->m_id;
532  if (v == w)
533  {
534  vertId[k] = tmpVertId[j];
535  break;
536  }
537  }
538  }
539  }
540  else if (el->GetConf().m_e == LibUtilities::ePrism)
541  {
542  boost::shared_ptr<Prism> pr =
543  boost::static_pointer_cast<Prism>(el);
544  if (pr->m_orientation == 1)
545  {
546  swap(vertId[2], vertId[1]);
547  swap(vertId[0], vertId[1]);
548  }
549  else if (pr->m_orientation == 2)
550  {
551  swap(vertId[0], vertId[1]);
552  swap(vertId[2], vertId[1]);
553  }
554  }
555 
556  HOSurfSharedPtr hs =
557  boost::shared_ptr<HOSurf>(new HOSurf(vertId));
558  // Find vertex combination in hoData.
559  hoIt = hoData[word].find(hs);
560 
561  if (hoIt == hoData[word].end())
562  {
563  cerr << "Unable to find high-order surface data "
564  << "for element id " << elId + 1 << endl;
565  abort();
566  }
567 
568  // Depending on order of vertices in rea file, surface
569  // information may need to be rotated or reflected.
570  HOSurfSharedPtr surf = *hoIt;
571  surf->Align(vertId);
572 
573  // Finally, add high order data to appropriate face.
574  int Ntot = (*hoIt)->surfVerts.size();
575  int N = ((int)sqrt(8.0 * Ntot + 1.0) - 1) / 2;
576  EdgeSharedPtr edge;
577 
578  // Apply high-order map to convert face data to Nektar++
579  // ordering (vertices->edges->internal).
580  vector<NodeSharedPtr> tmpVerts = (*hoIt)->surfVerts;
581  for (j = 0; j < tmpVerts.size(); ++j)
582  {
583  (*hoIt)->surfVerts[hoMap[j]] = tmpVerts[j];
584  }
585 
586  for (j = 0; j < tmpVerts.size(); ++j)
587  {
588  NodeSharedPtr a = (*hoIt)->surfVerts[j];
589  }
590 
591  vector<int> faceVertIds(3);
592  faceVertIds[0] = f->m_vertexList[0]->m_id;
593  faceVertIds[1] = f->m_vertexList[1]->m_id;
594  faceVertIds[2] = f->m_vertexList[2]->m_id;
595 
596  for (j = 0; j < f->m_edgeList.size(); ++j)
597  {
598  edge = f->m_edgeList[j];
599 
600  // Skip over edges which have already been
601  // populated.
602  if (edge->m_edgeNodes.size() > 0)
603  {
604  continue;
605  }
606 
607  // edge->m_edgeNodes.clear();
608  edge->m_curveType = LibUtilities::eGaussLobattoLegendre;
609 
610  for (int k = 0; k < N - 2; ++k)
611  {
612  edge->m_edgeNodes.push_back(
613  (*hoIt)->surfVerts[3 + j * (N - 2) + k]);
614  }
615 
616  // Nodal triangle data is always
617  // counter-clockwise. Add this check to reorder
618  // where necessary.
619  if (edge->m_n1->m_id != faceVertIds[j])
620  {
621  reverse(edge->m_edgeNodes.begin(),
622  edge->m_edgeNodes.end());
623  }
624  }
625 
626  // Insert interior face curvature.
627  f->m_curveType = LibUtilities::eNodalTriElec;
628  for (int j = 3 + 3 * (N - 2); j < Ntot; ++j)
629  {
630  f->m_faceNodes.push_back((*hoIt)->surfVerts[j]);
631  }
632  }
633  }
634  }
635 
636  // -- Process fluid boundary conditions.
637 
638  // Define a fairly horrendous map: key is the condition ID, the
639  // value is a vector of pairs of composites and element
640  // types. Essentially this map takes conditions -> composites for
641  // each element type.
642  map<int, vector<pair<int, LibUtilities::ShapeType> > > surfaceCompMap;
643 
644  // Skip boundary conditions line.
645  getline(m_mshFile, line);
646  getline(m_mshFile, line);
647 
648  int nSurfaces = 0;
649 
650  while (true)
651  {
652  getline(m_mshFile, line);
653 
654  // Break out of loop at end of boundary conditions section.
655  if (line.find("*") != string::npos || m_mshFile.eof() ||
656  line.length() == 0)
657  {
658  break;
659  }
660 
661  // Read boundary type, element ID and face ID.
662  char bcType;
663  int elId, faceId;
664  s.clear();
665  s.str(line);
666  s >> bcType >> elId >> faceId;
667  faceId--;
668  elId = elMap[elId - 1];
669 
670  vector<string> vals;
671  vector<ConditionType> type;
673 
674  ElementSharedPtr elm = m_mesh->m_element[m_mesh->m_spaceDim][elId];
675 
676  // Ignore BCs for undefined edges/faces
677  if ((elm->GetDim() == 2 && faceId >= elm->GetEdgeCount()) ||
678  (elm->GetDim() == 3 && faceId >= elm->GetFaceCount()))
679  {
680  continue;
681  }
682 
683  // First character on each line describes type of BC. Currently
684  // only support V, W, and O. In this switch statement we
685  // construct the quantities needed to search for the condition.
686  switch (bcType)
687  {
688  // Wall boundary.
689  case 'W':
690  {
691  if (scalar)
692  {
693  vals.push_back("0");
694  type.push_back(eDirichlet);
695  }
696  else
697  {
698  for (i = 0; i < m_mesh->m_fields.size() - 1; ++i)
699  {
700  vals.push_back("0");
701  type.push_back(eDirichlet);
702  }
703  // Set high-order boundary condition for wall.
704  vals.push_back("0");
705  type.push_back(eHOPCondition);
706  }
707  break;
708  }
709 
710  // Velocity boundary condition (either constant or dependent
711  // upon x,y,z).
712  case 'V':
713  case 'v':
714  {
715  if (scalar)
716  {
717  getline(m_mshFile, line);
718  size_t p = line.find_first_of('=');
719  vals.push_back(
720  boost::algorithm::trim_copy(line.substr(p + 1)));
721  type.push_back(eDirichlet);
722  }
723  else
724  {
725  for (i = 0; i < m_mesh->m_fields.size() - 1; ++i)
726  {
727  getline(m_mshFile, line);
728  size_t p = line.find_first_of('=');
729  vals.push_back(
730  boost::algorithm::trim_copy(line.substr(p + 1)));
731  type.push_back(eDirichlet);
732  }
733  // Set high-order boundary condition for Dirichlet
734  // condition.
735  vals.push_back("0");
736  type.push_back(eHOPCondition);
737  }
738  break;
739  }
740 
741  // Natural outflow condition (default value = 0.0?)
742  case 'O':
743  {
744  if (scalar)
745  {
746  vals.push_back("0");
747  type.push_back(eNeumann);
748  }
749  else
750  {
751  for (i = 0; i < m_mesh->m_fields.size(); ++i)
752  {
753  vals.push_back("0");
754  type.push_back(eNeumann);
755  }
756  // Set zero Dirichlet condition for outflow.
757  type[m_mesh->m_fields.size() - 1] = eDirichlet;
758  }
759  break;
760  }
761 
762  // Ignore unsupported BCs
763  case 'P':
764  case 'E':
765  case 'F':
766  case 'f':
767  case 'N':
768  continue;
769  break;
770 
771  default:
772  cerr << "Unknown boundary condition type " << line[0] << endl;
773  abort();
774  }
775 
776  // Populate condition information.
777  c->field = m_mesh->m_fields;
778  c->type = type;
779  c->value = vals;
780 
781  // Now attempt to find this boundary condition inside
782  // m_mesh->condition. This is currently a linear search and should
783  // probably be made faster!
785  bool found = false;
786  for (it = m_mesh->m_condition.begin(); it != m_mesh->m_condition.end();
787  ++it)
788  {
789  if (c == it->second)
790  {
791  found = true;
792  break;
793  }
794  }
795 
796  int compTag, conditionId;
797  ElementSharedPtr surfEl;
798 
799  // Create element for face (3D) or segment (2D). At the moment
800  // this is a bit of a hack since high-order nodes are not
801  // copied, so some output modules (e.g. Gmsh) will not output
802  // correctly.
803  if (elm->GetDim() == 3)
804  {
805  // 3D elements may have been reoriented, so face IDs will
806  // change.
807  if (elm->GetConf().m_e == LibUtilities::ePrism && faceId % 2 == 0)
808  {
809  boost::shared_ptr<Prism> p =
810  boost::dynamic_pointer_cast<Prism>(elm);
811  if (p->m_orientation == 1)
812  {
813  faceId = (faceId + 2) % 6;
814  }
815  else if (p->m_orientation == 2)
816  {
817  faceId = (faceId + 4) % 6;
818  }
819  }
820  else if (elm->GetConf().m_e == LibUtilities::eTetrahedron)
821  {
822  boost::shared_ptr<Tetrahedron> t =
823  boost::dynamic_pointer_cast<Tetrahedron>(elm);
824  faceId = t->m_orientationMap[faceId];
825  }
826 
827  FaceSharedPtr f = elm->GetFace(faceId);
828  bool tri = f->m_vertexList.size() == 3;
829 
830  vector<NodeSharedPtr> nodeList;
831  nodeList.insert(nodeList.begin(),
832  f->m_vertexList.begin(),
833  f->m_vertexList.end());
834 
835  vector<int> tags;
836 
839  ElmtConfig conf(
840  seg, 1, true, true, false, LibUtilities::eGaussLobattoLegendre);
841  surfEl =
842  GetElementFactory().CreateInstance(seg, conf, nodeList, tags);
843 
844  // Copy high-order surface information from edges.
845  for (int i = 0; i < f->m_vertexList.size(); ++i)
846  {
847  surfEl->GetEdge(i)->m_edgeNodes = f->m_edgeList[i]->m_edgeNodes;
848  surfEl->GetEdge(i)->m_curveType = f->m_edgeList[i]->m_curveType;
849  }
850  }
851  else if (faceId < elm->GetEdgeCount())
852  {
853  EdgeSharedPtr f = elm->GetEdge(faceId);
854 
855  vector<NodeSharedPtr> nodeList;
856  nodeList.push_back(f->m_n1);
857  nodeList.push_back(f->m_n2);
858 
859  vector<int> tags;
860 
862  1,
863  true,
864  true,
865  false,
868  LibUtilities::eSegment, conf, nodeList, tags);
869  }
870 
871  if (!surfEl)
872  {
873  continue;
874  }
875 
876  LibUtilities::ShapeType surfElType = surfEl->GetConf().m_e;
877 
878  if (!found)
879  {
880  // If condition does not already exist, add to condition
881  // list, create new composite tag and put inside
882  // surfaceCompMap.
883  conditionId = m_mesh->m_condition.size();
884  compTag = nComposite;
885  c->m_composite.push_back(compTag);
886  m_mesh->m_condition[conditionId] = c;
887 
888  surfaceCompMap[conditionId].push_back(
889  pair<int, LibUtilities::ShapeType>(nComposite, surfElType));
890 
891  nComposite++;
892  }
893  else
894  {
895  // Otherwise find existing composite inside surfaceCompMap.
896  map<int, vector<pair<int, LibUtilities::ShapeType> > >::iterator
897  it2;
898  it2 = surfaceCompMap.find(it->first);
899 
900  found = false;
901  if (it2 == surfaceCompMap.end())
902  {
903  // This should never happen!
904  cerr << "Unable to find condition!" << endl;
905  abort();
906  }
907 
908  for (j = 0; j < it2->second.size(); ++j)
909  {
910  pair<int, LibUtilities::ShapeType> tmp = it2->second[j];
911  if (tmp.second == surfElType)
912  {
913  found = true;
914  compTag = tmp.first;
915  break;
916  }
917  }
918 
919  // If no pairs were found, then this condition contains
920  // multiple element types (i.e. both triangles and
921  // quads). Create another composite for the new shape type
922  // and insert into the map.
923  if (!found)
924  {
925  it2->second.push_back(
926  pair<int, LibUtilities::ShapeType>(nComposite, surfElType));
927  compTag = nComposite;
928  m_mesh->m_condition[it->first]->m_composite.push_back(compTag);
929  nComposite++;
930  }
931 
932  conditionId = it->first;
933  }
934 
935  // Insert composite tag into element and insert element into
936  // mesh.
937  vector<int> existingTags = surfEl->GetTagList();
938 
939  existingTags.insert(existingTags.begin(), compTag);
940  surfEl->SetTagList(existingTags);
941  surfEl->SetId(nSurfaces);
942 
943  m_mesh->m_element[surfEl->GetDim()].push_back(surfEl);
944  nSurfaces++;
945  }
946 
947  m_mshFile.reset();
948 
949  // -- Process rest of mesh.
950  ProcessEdges();
951  ProcessFaces();
952  ProcessElements();
954 }
A 3-dimensional four-faced element.
Definition: Tetrahedron.h:50
Basic information about an element.
Definition: ElementConfig.h:50
io::filtering_istream m_mshFile
Input stream.
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
A 3-dimensional five-faced element (2 triangles, 3 quadrilaterals).
Definition: Prism.h:50
boost::shared_ptr< HOSurf > HOSurfSharedPtr
Definition: HOAlignment.h:181
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
std::map< int, int > hoMap
Definition: InputNek.h:92
NEKMESHUTILS_EXPORT void OpenStream()
Open a file for input.
std::map< std::string, NekMeshUtils::HOSurfSet > hoData
Definition: InputNek.h:87
NEKMESHUTILS_EXPORT NodeSharedPtr GetVertex(unsigned int i) const
Access a vertex node.
Definition: Element.h:100
virtual NEKMESHUTILS_EXPORT void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
virtual NEKMESHUTILS_EXPORT void ProcessElements()
Generate element IDs.
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
std::map< std::string, ConfigOption > m_config
List of configuration values.
boost::shared_ptr< Condition > ConditionSharedPtr
Definition: Mesh.h:82
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
Definition: Edge.h:136
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
int GetNnodes(LibUtilities::ShapeType elType)
Definition: InputNek.cpp:1124
std::map< std::string, std::pair< NekCurve, std::string > > curveTags
Definition: InputNek.h:82
virtual NEKMESHUTILS_EXPORT void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
boost::shared_ptr< Face > FaceSharedPtr
Definition: Face.h:153
HOTriangle< NodeSharedPtr > HOSurf
Definition: HOAlignment.h:180
virtual NEKMESHUTILS_EXPORT void ProcessComposites()
Generate composites.
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:52
2D Nodal Electrostatic Points on a Triangle
Definition: PointsType.h:70

Member Data Documentation

ModuleKey Nektar::Utilities::InputNek::className
static
Initial value:

ModuleKey for class.

Definition at line 73 of file InputNek.h.

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

Referenced by LoadHOSurfaces(), and Process().

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

Referenced by LoadHOSurfaces(), and Process().