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. More...
 
- 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. More...
 
- Public Member Functions inherited from Nektar::Utilities::Module
 Module (FieldSharedPtr p_f)
 
virtual void Process (po::variables_map &vm)=0
 
virtual std::string GetModuleName ()=0
 
void RegisterConfig (string key, string value)
 Register a configuration option with a module. More...
 
void PrintConfig ()
 Print out all configuration options for a module. More...
 
void SetDefaults ()
 Sets default configuration options for those which have not been set. More...
 
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 (std::string key, std::string value)
 
void PrintConfig ()
 
void SetDefaults ()
 
MeshSharedPtr GetMesh ()
 
virtual void ProcessVertices ()
 Extract element vertices. More...
 
virtual void ProcessEdges (bool ReprocessEdges=true)
 Extract element edges. More...
 
virtual void ProcessFaces (bool ReprocessFaces=true)
 Extract element faces. More...
 
virtual void ProcessElements ()
 Generate element IDs. More...
 
virtual void ProcessComposites ()
 Generate composites. More...
 
virtual void ClearElementLinks ()
 

Static Public Member Functions

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

Static Public Attributes

static 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, HOSurfSethoData
 
std::map< int, int > hoMap
 

Additional Inherited Members

- Protected Member Functions inherited from Nektar::Utilities::InputModule
void PrintSummary ()
 Print summary of elements. More...
 
void PrintSummary ()
 Print summary of elements. More...
 
- Protected Member Functions inherited from Nektar::Utilities::Module
 Module ()
 
void ReorderPrisms (PerMap &perFaces)
 Reorder node IDs so that prisms and tetrahedra are aligned correctly. More...
 
void PrismLines (int prism, PerMap &perFaces, std::set< int > &prismsDone, std::vector< ElementSharedPtr > &line)
 
- Protected Attributes inherited from Nektar::Utilities::InputModule
set< string > m_allowedFiles
 
std::ifstream m_mshFile
 Input stream. More...
 
- Protected Attributes inherited from Nektar::Utilities::Module
FieldSharedPtr m_f
 Field object. More...
 
map< string, ConfigOptionm_config
 List of configuration values. More...
 
bool m_requireEquiSpaced
 
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 ( MeshSharedPtr  p_m)

Definition at line 60 of file InputNek.cpp.

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

Definition at line 64 of file InputNek.cpp.

65 {
66 }

Member Function Documentation

static ModuleSharedPtr Nektar::Utilities::InputNek::create ( 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 1076 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().

1077 {
1078  int nNodes = 0;
1079 
1080  switch (InputNekEntity)
1081  {
1082  case LibUtilities::ePoint:
1083  nNodes = 1;
1084  break;
1086  nNodes = 2;
1087  break;
1089  nNodes = 3;
1090  break;
1092  nNodes = 4;
1093  break;
1095  nNodes = 4;
1096  break;
1098  nNodes = 5;
1099  break;
1100  case LibUtilities::ePrism:
1101  nNodes = 6;
1102  break;
1104  nNodes = 8;
1105  break;
1106  default:
1107  cerr << "unknown Nektar element type" << endl;
1108  }
1109 
1110  return nNodes;
1111 }
void Nektar::Utilities::InputNek::LoadHOSurfaces ( )
private

Load high order surface information from hsf file.

Definition at line 911 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().

912 {
913  map<string, pair<NekCurve, string> >::iterator it;
914  int nodeId = m_mesh->GetNumEntities();
915 
916  for (it = curveTags.begin(); it != curveTags.end(); ++it)
917  {
918  ifstream hsf;
919  string line, fileName = it->second.second;
920  size_t pos;
921  int N, Nface, dot;
922 
923  if (it->second.first != eFile)
924  {
925  continue;
926  }
927 
928  // Replace fro extension with hsf.
929  dot = fileName.find_last_of('.');
930  fileName = fileName.substr(0, dot);
931  fileName += ".hsf";
932 
933  // Open hsf file.
934  hsf.open(fileName.c_str());
935  if (!hsf.is_open())
936  {
937  cerr << "Could not open surface file " << fileName << endl;
938  abort();
939  }
940 
941  // Read in header line; determine element order, number of faces
942  // from this line.
943  getline(hsf, line);
944  pos = line.find("=");
945  if (pos == string::npos)
946  {
947  cerr << "hsf header error: cannot read number of "
948  << "nodal points." << endl;
949  abort();
950  }
951  line = line.substr(pos + 1);
952  stringstream ss(line);
953  ss >> N;
954 
955  pos = line.find("=");
956  if (pos == string::npos)
957  {
958  cerr << "hsf header error: cannot read number of "
959  << "faces." << endl;
960  abort();
961  }
962  line = line.substr(pos + 1);
963  ss.clear();
964  ss.str(line);
965  ss >> Nface;
966 
967  int Ntot = N * (N + 1) / 2;
968 
969  // Skip a line, then read in r,s positions inside the next
970  // comments.
971  Array<OneD, NekDouble> r(Ntot), s(Ntot);
972  getline(hsf, line);
973 
974  for (int i = 0; i < 2; ++i)
975  {
976  string word;
977 
978  getline(hsf, line);
979  ss.clear();
980  ss.str(line);
981  ss >> word;
982 
983  if (word != "#")
984  {
985  cerr << "hsf header error: cannot read in "
986  << "r/s points" << endl;
987  abort();
988  }
989 
990  for (int j = 0; j < Ntot; ++j)
991  {
992  ss >> (i == 0 ? r[j] : s[j]);
993  }
994  }
995 
996  // Generate electrostatic points so that re-mapping array can
997  // be constructed.
998  Array<OneD, NekDouble> rp(Ntot), sp(Ntot);
999  LibUtilities::PointsKey elec(N, LibUtilities::eNodalTriElec);
1000  LibUtilities::PointsManager()[elec]->GetPoints(rp, sp);
1001 
1002  // Expensively construct remapping array nodemap. This will
1003  // map nodal ordering from hsf order to Nektar++ ordering
1004  // (i.e. vertices followed by edges followed by interior
1005  // points.)
1006  for (int i = 0; i < Ntot; ++i)
1007  {
1008  for (int j = 0; j < Ntot; ++j)
1009  {
1010  if (fabs(r[i] - rp[j]) < 1e-5 && fabs(s[i] - sp[j]) < 1e-5)
1011  {
1012  hoMap[i] = j;
1013  break;
1014  }
1015  }
1016  }
1017 
1018  // Skip variables line
1019  getline(hsf, line);
1020 
1021  // Read in nodal points for each face.
1022  map<int, vector<NodeSharedPtr> > faceMap;
1023  for (int i = 0; i < Nface; ++i)
1024  {
1025  getline(hsf, line);
1026  vector<NodeSharedPtr> faceNodes(Ntot);
1027  for (int j = 0; j < Ntot; ++j, ++nodeId)
1028  {
1029  double x, y, z;
1030  getline(hsf, line);
1031  ss.clear();
1032  ss.str(line);
1033  ss >> x >> y >> z;
1034  faceNodes[j] = NodeSharedPtr(new Node(nodeId, x, y, z));
1035  }
1036  // Skip over tecplot connectivity information.
1037  for (int j = 0; j < (N - 1) * (N - 1); ++j)
1038  {
1039  getline(hsf, line);
1040  }
1041  faceMap[i] = faceNodes;
1042  }
1043 
1044  // Finally, read in connectivity information to set up after
1045  // reading rea file.
1046  getline(hsf, line);
1047  for (int i = 0; i < Nface; ++i)
1048  {
1049  string tmp;
1050  int fid;
1051  vector<int> nodeIds(3);
1052 
1053  getline(hsf, line);
1054  ss.clear();
1055  ss.str(line);
1056  ss >> tmp >> fid >> nodeIds[0] >> nodeIds[1] >> nodeIds[2];
1057 
1058  if (tmp != "#")
1059  {
1060  cerr << "Unable to read hsf connectivity information." << endl;
1061  abort();
1062  }
1063 
1064  hoData[it->first].insert(
1065  HOSurfSharedPtr(new HOSurf(nodeIds, faceMap[i])));
1066  }
1067 
1068  hsf.close();
1069  }
1070 }
boost::shared_ptr< HOSurf > HOSurfSharedPtr
Definition: Triangle.h:216
MeshSharedPtr m_mesh
Mesh object.
std::map< std::string, HOSurfSet > hoData
Definition: InputNek.h:87
Represents a point in the domain.
Definition: Node.h:60
std::map< int, int > hoMap
Definition: InputNek.h:92
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: Triangle.h:215
2D Nodal Electrostatic Points on a Triangle
Definition: PointsType.h:68
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 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::Utilities::Module::m_mesh, Nektar::Utilities::InputModule::m_mshFile, Nektar::NekMeshUtils::Tetrahedron::m_orientationMap, Nektar::Utilities::InputModule::OpenStream(), Nektar::Utilities::Module::ProcessComposites(), Nektar::Utilities::Module::ProcessEdges(), Nektar::Utilities::Module::ProcessElements(), Nektar::Utilities::Module::ProcessFaces(), and Nektar::NekMeshUtils::surf.

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

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, 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().