Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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  m_config["scalar"] = ConfigOption(
63  true, "0", "If defined then assume input rea is for scalar problem");
64 }
65 
67 {
68 }
69 
70 /**
71  * @brief Processes Nektar file format.
72  *
73  * Nektar sessions are defined by rea files, and contain sections defining a DNS
74  * simulation in a specific order. The converter only reads mesh information,
75  * curve information if it exists and boundary information.
76  *
77  * @param pFilename Filename of Nektar session file to read.
78  */
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 }
955 
956 /**
957  * Load high order surface information from hsf file.
958  */
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);
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 }
1119 
1120 /**
1121  * This routine aids the reading of Nektar files only; it returns the
1122  * number of nodes for a given entity typw.
1123  */
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 }
1160 }
1161 }
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
STL namespace.
pair< ModuleType, string > ModuleKey
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.
Represents a command-line configuration option.
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
PointsManagerT & PointsManager(void)
Defines a specification for a set of points.
Definition: Points.h:58
std::map< std::string, ConfigOption > m_config
List of configuration values.
boost::shared_ptr< Condition > ConditionSharedPtr
Definition: Mesh.h:82
virtual void Process()
Processes Nektar file format.
Definition: InputNek.cpp:79
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
Definition: Edge.h:136
Abstract base class for input modules.
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:147
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
std::pair< ModuleType, std::string > ModuleKey
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
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215