Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
InputNek.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: InputNek.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Nektar file format converter.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <map>
37 #include <vector>
38 #include <sstream>
39 #include <string>
40 
41 #include <boost/algorithm/string.hpp>
46 
47 #include "InputNek.h"
48 
49 using namespace std;
50 using namespace Nektar::NekMeshUtils;
51 
52 namespace Nektar
53 {
54 namespace Utilities
55 {
56 
57 ModuleKey InputNek::className = GetModuleFactory().RegisterCreatorFunction(
58  ModuleKey(eInputModule, "rea"), InputNek::create, "Reads Nektar rea file.");
59 
60 InputNek::InputNek(MeshSharedPtr m) : InputModule(m)
61 {
62 }
63 
65 {
66 }
67 
68 /**
69  * @brief Processes Nektar file format.
70  *
71  * Nektar sessions are defined by rea files, and contain sections defining a DNS
72  * simulation in a specific order. The converter only reads mesh information,
73  * curve information if it exists and boundary information.
74  *
75  * @param pFilename Filename of Nektar session file to read.
76  */
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 }
907 
908 /**
909  * Load high order surface information from hsf file.
910  */
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);
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 }
1071 
1072 /**
1073  * This routine aids the reading of Nektar files only; it returns the
1074  * number of nodes for a given entity typw.
1075  */
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 }
1112 }
1113 }
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.
pair< ModuleType, string > ModuleKey
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
STL namespace.
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
PointsManagerT & PointsManager(void)
Defines a specification for a set of points.
Definition: Points.h:58
boost::shared_ptr< Condition > ConditionSharedPtr
Definition: Mesh.h:78
virtual void ProcessComposites()
Generate composites.
virtual void Process()
Processes Nektar file format.
Definition: InputNek.cpp:77
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
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Definition: Mesh.h:137
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
ModuleFactory & GetModuleFactory()
Abstract base class for input modules.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215