Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
InputStarTec.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: InputStarTec.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: Tecplot file converter.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <boost/algorithm/string.hpp>
37 
40 
41 #include "InputStarTec.h"
42 
43 using namespace std;
44 using namespace Nektar::NekMeshUtils;
45 
46 namespace Nektar
47 {
48 namespace Utilities
49 {
50 
51 ModuleKey InputTec::className = GetModuleFactory().RegisterCreatorFunction(
52  ModuleKey(eInputModule, "dat"),
53  InputTec::create,
54  "Reads Tecplot polyhedron ascii format converted from Star CCM (.dat).");
55 
56 InputTec::InputTec(MeshSharedPtr m) : InputModule(m)
57 {
58 }
59 
61 {
62 }
63 
64 /**
65  * Tecplot file Polyhedron format contains a list of nodes, a
66  * node count per face, the node ids, Element ids that are on
67  * the left of each face and Element ids which are on the
68  * right of each face. There are then a series of zone of each
69  * surface. In the case of a surface the number of nodes is
70  * not provided indicating it is a 2D zone
71  *
72  * @param pFilename Filename of Tecplot file to read.
73  */
75 {
76  m_mesh->m_expDim = 3;
77  m_mesh->m_spaceDim = 3;
78 
79  if (m_mesh->m_verbose)
80  {
81  cout << "InputStarTec: Start reading file..." << endl;
82  }
83 
84  string line, word;
85 
86  // Open the file stream.
87  OpenStream();
88 
89  int nComposite = 0;
90 
91  // read first zone (Hopefully 3D)
92  while (!m_mshFile.eof())
93  {
94  getline(m_mshFile, line);
95  if (line.find("ZONE") != string::npos)
96  {
97  ReadZone(nComposite);
98  break;
99  }
100  }
101 
102  // read remaining 2D zones
103  while (!m_mshFile.eof())
104  {
105  if (line.find("ZONE") != string::npos)
106  {
107  ReadZone(nComposite);
108  }
109  }
110 
111  PrintSummary();
112  m_mshFile.close();
113 
114  ProcessEdges();
115  ProcessFaces();
116  ProcessElements();
118 }
119 
120 void InputTec::ReadZone(int &nComposite)
121 {
122  int i;
123  string line, tag;
124  int nfaces, nnodes, nelements;
125  int start, end;
126  stringstream s;
127  NekDouble value;
128  streampos pos;
129  static int zcnt = 1;
130 
131  // Read Zone Header
132  nnodes = nfaces = nelements = 0;
133  while (!m_mshFile.eof())
134  {
135  pos = m_mshFile.tellg();
136 
137  getline(m_mshFile, line);
138 
139  boost::to_upper(line);
140 
141  // cehck to see if readable data.
142  if (sscanf(line.c_str(), "%lf", &value) == 1)
143  {
144  m_mshFile.seekg(pos);
145  break;
146  }
147 
148  if ((line.find("NODES") != string::npos) &&
149  (line.find("TOTALNUMFACENODES") == string::npos))
150  {
151  s.clear();
152  s.str(line);
153 
154  tag = s.str();
155  start = tag.find("NODES=");
156  end = tag.find_first_of(',', start);
157  nnodes = atoi(tag.substr(start + 6, end).c_str());
158  }
159 
160  if ((line.find("FACES") != string::npos) &&
161  (line.find("NUMCONNECTEDBOUNDARYFACES") == string::npos))
162  {
163  s.clear();
164  s.str(line);
165 
166  tag = s.str();
167  start = tag.find("FACES=");
168  end = tag.find_first_of(',', start);
169  nfaces = atoi(tag.substr(start + 6, end).c_str());
170  }
171 
172  if (line.find("ELEMENTS") != string::npos)
173  {
174  s.clear();
175  s.str(line);
176 
177  tag = s.str();
178  start = tag.find("ELEMENTS=");
179  end = tag.find_first_of(',', start);
180  nelements = atoi(tag.substr(start + 9, end).c_str());
181  }
182 
183  if (line.find("ZONETYPE") != string::npos)
184  {
185  s.clear();
186  s.str(line);
187 
188  if ((line.find("FEPOLYGON") == string::npos) &&
189  (line.find("FEPOLYHEDRON") == string::npos))
190  {
191  ASSERTL0(false,
192  "Routine only set up for FEPolygon or FEPolyhedron");
193  }
194  }
195  }
196  if (!nnodes) // No zone found
197  {
198  return;
199  }
200 
201  cout << "Setting up zone " << zcnt++;
202 
203  vector<NekDouble> x, y, z;
204 
205  // Read in Nodes
206  for (i = 0; i < nnodes; ++i)
207  {
208  m_mshFile >> value;
209  x.push_back(value);
210  }
211 
212  for (i = 0; i < nnodes; ++i)
213  {
214  m_mshFile >> value;
215  y.push_back(value);
216  }
217 
218  for (i = 0; i < nnodes; ++i)
219  {
220  m_mshFile >> value;
221  z.push_back(value);
222  }
223 
224  std::vector<NodeSharedPtr> Nodes;
225  for (i = 0; i < nnodes; ++i)
226  {
227  Nodes.push_back(boost::shared_ptr<Node>(new Node(i, x[i], y[i], z[i])));
228  }
229 
230  // Read Node count per face
231  getline(m_mshFile, line);
232  if (line.find("node count per face") == string::npos)
233  {
234  if (line.find("face nodes") == string::npos)
235  {
236  getline(m_mshFile, line);
237  }
238  }
239 
240  s.clear();
241  s.str(line);
242 
243  vector<int> Nodes_per_face;
244  if (line.find("node count per face") != string::npos)
245  {
246  int nodes;
247  for (i = 0; i < nfaces; ++i)
248  {
249  m_mshFile >> nodes;
250  ASSERTL0(nodes <= 4,
251  "Can only handle meshes with "
252  "up to four nodes per face");
253  Nodes_per_face.push_back(nodes);
254  }
255  // Read next line
256  getline(m_mshFile, line);
257  }
258 
259  // Read face nodes;
260  if (line.find("face nodes") == string::npos)
261  {
262  getline(m_mshFile, line);
263  }
264  s.clear();
265  s.str(line);
266 
267  vector<vector<int> > FaceNodes;
268 
269  if (line.find("face nodes") != string::npos)
270  {
271 
272  for (i = 0; i < nfaces; ++i)
273  {
274  // check to see if Nodes_per_face is defined and
275  // if not assume 2 nodes for 2D case
276  int nodes = (Nodes_per_face.size()) ? Nodes_per_face[i] : 2;
277 
278  int nodeID;
279  vector<int> Fnodes;
280  for (int j = 0; j < nodes; ++j)
281  {
282 
283  m_mshFile >> nodeID;
284 
285  Fnodes.push_back(nodeID - 1);
286  }
287 
288  FaceNodes.push_back(Fnodes);
289  }
290  }
291  else
292  {
293  ASSERTL0(false, "Failed to find face node section");
294  }
295 
296  // Read left elements
297  Array<OneD, vector<int> > ElementFaces(nelements);
298 
299  // check to see if next line contains left elements
300  getline(m_mshFile, line);
301  if (line.find("left elements") == string::npos)
302  {
303  getline(m_mshFile, line);
304  }
305 
306  if (line.find("left elements") != string::npos)
307  {
308  int elmtID;
309 
310  for (i = 0; i < nfaces; ++i)
311  {
312  m_mshFile >> elmtID;
313 
314  if (elmtID > 0)
315  {
316  ElementFaces[elmtID - 1].push_back(i);
317  }
318  }
319  }
320  else
321  {
322  ASSERTL0(false, "Left element not found");
323  }
324 
325  // check to see if next line contains right elements
326  getline(m_mshFile, line);
327  if (line.find("right elements") == string::npos)
328  {
329  getline(m_mshFile, line);
330  }
331 
332  if (line.find("right elements") != string::npos)
333 
334  {
335  int elmtID;
336 
337  for (i = 0; i < nfaces; ++i)
338  {
339  m_mshFile >> elmtID;
340 
341  if (elmtID > 0)
342  {
343  ElementFaces[elmtID - 1].push_back(i);
344  }
345  }
346 
347  // read to end of line
348  getline(m_mshFile, line);
349  }
350  else
351  {
352  ASSERTL0(false, "Left element not found");
353  }
354 
355  if (Nodes_per_face.size()) // 3D Zone
356  {
357  cout << " (3D) " << endl;
358 
359  // Reset node ordering so that all prism faces have
360  // consistent numbering for singular vertex re-ordering
361  ResetNodes(Nodes, ElementFaces, FaceNodes);
362 
363  m_mesh->m_node = Nodes;
364 
365  // create Prisms/Pyramids first
366  for (i = 0; i < nelements; ++i)
367  {
368  if (ElementFaces[i].size() > 4)
369  {
370  GenElement3D(
371  Nodes, i, ElementFaces[i], FaceNodes, nComposite, true);
372  }
373  }
374 
375  nComposite++;
376 
377  // create Tets second
378  for (i = 0; i < nelements; ++i)
379  {
380  if (ElementFaces[i].size() == 4)
381  {
382  GenElement3D(
383  Nodes, i, ElementFaces[i], FaceNodes, nComposite, true);
384  }
385  }
386  nComposite++;
387 
388  ProcessVertices();
389  }
390  else // 2D Zone
391  {
392  cout << " (2D)" << endl;
393 
394  // find ids of VertNodes from m_mesh->m_vertexSet so that we can
395  // identify
396  for (i = 0; i < Nodes.size(); ++i)
397  {
398  NodeSet::iterator it = m_mesh->m_vertexSet.find(Nodes[i]);
399 
400  if (it == m_mesh->m_vertexSet.end())
401  {
402  ASSERTL0(false, "Failed to find face vertex in 3D list");
403  }
404  else
405  {
406  Nodes[i] = *it;
407  }
408  }
409 
410  for (i = 0; i < nelements; ++i)
411  {
412  GenElement2D(Nodes, i, ElementFaces[i], FaceNodes, nComposite);
413  }
414 
415  nComposite++;
416  }
417 }
418 
419 static void PrismLineFaces(int prismid,
420  map<int, int> &facelist,
421  vector<vector<int> > &FacesToPrisms,
422  vector<vector<int> > &PrismsToFaces,
423  vector<bool> &PrismDone);
424 
425 void InputTec::ResetNodes(vector<NodeSharedPtr> &Vnodes,
426  Array<OneD, vector<int> > &ElementFaces,
427  vector<vector<int> > &FaceNodes)
428 {
429  int i, j;
430  Array<OneD, int> NodeReordering(Vnodes.size(), -1);
431  int face1_map[3] = {0, 1, 4};
432  int face3_map[3] = {3, 2, 5};
433  int nodeid = 0;
434  map<int, bool> FacesRenumbered;
435 
436  // Determine Prism triangular face connectivity.
437  vector<vector<int> > FaceToPrisms(FaceNodes.size());
438  vector<vector<int> > PrismToFaces(ElementFaces.num_elements());
439  map<int, int> Prisms;
440  map<int, int>::iterator PrismIt;
441 
442  // generate map of prism-faces to prisms and prism to
443  // triangular-faces as well as ids of each prism.
444  for (i = 0; i < ElementFaces.num_elements(); ++i)
445  {
446  // Find Prism (and pyramids!).
447  if (ElementFaces[i].size() == 5)
448  {
449  vector<int> LocTriFaces;
450  // Find triangular faces
451  for (j = 0; j < ElementFaces[i].size(); ++j)
452  {
453  if (FaceNodes[ElementFaces[i][j]].size() == 3)
454  {
455  LocTriFaces.push_back(j);
456  }
457  }
458 
459  if (LocTriFaces.size() == 2) // prism otherwise a pyramid
460  {
461  Prisms[i] = i;
462 
463  PrismToFaces[i].push_back(ElementFaces[i][LocTriFaces[0]]);
464  PrismToFaces[i].push_back(ElementFaces[i][LocTriFaces[1]]);
465 
466  FaceToPrisms[ElementFaces[i][LocTriFaces[0]]].push_back(i);
467  FaceToPrisms[ElementFaces[i][LocTriFaces[1]]].push_back(i);
468  }
469  }
470  }
471 
472  vector<bool> FacesDone(FaceNodes.size(), false);
473  vector<bool> PrismDone(ElementFaces.num_elements(), false);
474 
475  // For every prism find the list of prismatic elements
476  // that represent an aligned block of cells. Then renumber
477  // these blocks consecutativiesly
478  for (PrismIt = Prisms.begin(); PrismIt != Prisms.end(); ++PrismIt)
479  {
480  int elmtid = PrismIt->first;
481  map<int, int> facelist;
483 
484  if (PrismDone[elmtid])
485  {
486  continue;
487  }
488  else
489  {
490  // Generate list of faces in list
492  elmtid, facelist, FaceToPrisms, PrismToFaces, PrismDone);
493 
494  // loop over faces and number vertices of associated prisms.
495  for (faceIt = facelist.begin(); faceIt != facelist.end(); faceIt++)
496  {
497  int faceid = faceIt->second;
498 
499  for (i = 0; i < FaceToPrisms[faceid].size(); ++i)
500  {
501  int prismid = FaceToPrisms[faceid][i];
502 
503  if ((FacesDone[PrismToFaces[prismid][0]] == true) &&
504  (FacesDone[PrismToFaces[prismid][1]] == true))
505  {
506  continue;
507  }
508 
509  Array<OneD, int> Nodes =
510  SortFaceNodes(Vnodes, ElementFaces[prismid], FaceNodes);
511 
512  if ((FacesDone[PrismToFaces[prismid][0]] == false) &&
513  (FacesDone[PrismToFaces[prismid][1]] == false))
514  {
515  // number all nodes consecutive since
516  // already correctly re-arranged.
517  for (i = 0; i < 3; ++i)
518  {
519  if (NodeReordering[Nodes[face1_map[i]]] == -1)
520  {
521  NodeReordering[Nodes[face1_map[i]]] = nodeid++;
522  }
523  }
524 
525  for (i = 0; i < 3; ++i)
526  {
527  if (NodeReordering[Nodes[face3_map[i]]] == -1)
528  {
529  NodeReordering[Nodes[face3_map[i]]] = nodeid++;
530  }
531  }
532  }
533  else if ((FacesDone[PrismToFaces[prismid][0]] == false) &&
534  (FacesDone[PrismToFaces[prismid][1]] == true))
535  {
536  // find node of highest id
537  int max_id1, max_id2;
538 
539  max_id1 = (NodeReordering[Nodes[face3_map[0]]] <
540  NodeReordering[Nodes[face3_map[1]]])
541  ? 1
542  : 0;
543  max_id2 = (NodeReordering[Nodes[face3_map[max_id1]]] <
544  NodeReordering[Nodes[face3_map[2]]])
545  ? 2
546  : max_id1;
547 
548  // add numbering according to order of
549  int id0 = (max_id1 == 1) ? 0 : 1;
550 
551  if (NodeReordering[Nodes[face1_map[id0]]] == -1)
552  {
553  NodeReordering[Nodes[face1_map[id0]]] = nodeid++;
554  }
555 
556  if (NodeReordering[Nodes[face1_map[max_id1]]] == -1)
557  {
558  NodeReordering[Nodes[face1_map[max_id1]]] =
559  nodeid++;
560  }
561 
562  if (NodeReordering[Nodes[face1_map[max_id2]]] == -1)
563  {
564  NodeReordering[Nodes[face1_map[max_id2]]] =
565  nodeid++;
566  }
567  }
568  else if ((FacesDone[PrismToFaces[prismid][0]] == true) &&
569  (FacesDone[PrismToFaces[prismid][1]] == false))
570  {
571  // find node of highest id
572  int max_id1, max_id2;
573 
574  max_id1 = (NodeReordering[Nodes[face1_map[0]]] <
575  NodeReordering[Nodes[face1_map[1]]])
576  ? 1
577  : 0;
578  max_id2 = (NodeReordering[Nodes[face1_map[max_id1]]] <
579  NodeReordering[Nodes[face1_map[2]]])
580  ? 2
581  : max_id1;
582 
583  // add numbering according to order of
584  int id0 = (max_id1 == 1) ? 0 : 1;
585 
586  if (NodeReordering[Nodes[face3_map[id0]]] == -1)
587  {
588  NodeReordering[Nodes[face3_map[id0]]] = nodeid++;
589  }
590 
591  if (NodeReordering[Nodes[face3_map[max_id1]]] == -1)
592  {
593  NodeReordering[Nodes[face3_map[max_id1]]] =
594  nodeid++;
595  }
596 
597  if (NodeReordering[Nodes[face3_map[max_id2]]] == -1)
598  {
599  NodeReordering[Nodes[face3_map[max_id2]]] =
600  nodeid++;
601  }
602  }
603  }
604  }
605  }
606  }
607 
608  // fill in any unset nodes at from other shapes
609  for (i = 0; i < NodeReordering.num_elements(); ++i)
610  {
611  if (NodeReordering[i] == -1)
612  {
613  NodeReordering[i] = nodeid++;
614  }
615  }
616 
617  ASSERTL1(nodeid == NodeReordering.num_elements(),
618  "Have not renumbered all nodes");
619 
620  // Renumbering successfull so resort nodes and faceNodes;
621  for (i = 0; i < FaceNodes.size(); ++i)
622  {
623  for (j = 0; j < FaceNodes[i].size(); ++j)
624  {
625  FaceNodes[i][j] = NodeReordering[FaceNodes[i][j]];
626  }
627  }
628 
629  vector<NodeSharedPtr> save(Vnodes);
630  for (i = 0; i < Vnodes.size(); ++i)
631  {
632  Vnodes[NodeReordering[i]] = save[i];
633  Vnodes[NodeReordering[i]]->SetID(NodeReordering[i]);
634  }
635 }
636 
637 static void PrismLineFaces(int prismid,
638  map<int, int> &facelist,
639  vector<vector<int> > &FaceToPrisms,
640  vector<vector<int> > &PrismToFaces,
641  vector<bool> &PrismDone)
642 {
643  if (PrismDone[prismid] == false)
644  {
645  PrismDone[prismid] = true;
646 
647  // Add faces0
648  int face = PrismToFaces[prismid][0];
649  facelist[face] = face;
650  for (int i = 0; i < FaceToPrisms[face].size(); ++i)
651  {
652  PrismLineFaces(FaceToPrisms[face][i],
653  facelist,
654  FaceToPrisms,
655  PrismToFaces,
656  PrismDone);
657  }
658 
659  // Add faces1
660  face = PrismToFaces[prismid][1];
661  facelist[face] = face;
662  for (int i = 0; i < FaceToPrisms[face].size(); ++i)
663  {
664  PrismLineFaces(FaceToPrisms[face][i],
665  facelist,
666  FaceToPrisms,
667  PrismToFaces,
668  PrismDone);
669  }
670  }
671 }
672 
673 void InputTec::GenElement2D(vector<NodeSharedPtr> &VertNodes,
674  int i,
675  vector<int> &ElementFaces,
676  vector<vector<int> > &FaceNodes,
677  int nComposite)
678 {
680  // set up Node list
681 
682  if (ElementFaces.size() == 3)
683  {
684  elType = LibUtilities::eTriangle;
685  }
686  else if (ElementFaces.size() == 4)
687  {
689  }
690  else
691  {
692  ASSERTL0(false, "Not set up for elements which are not Tets or Prism");
693  }
694 
695  // Create element tags
696  vector<int> tags;
697  tags.push_back(nComposite);
698 
699  // make unique node list
700  vector<NodeSharedPtr> nodeList;
701  Array<OneD, int> Nodes = SortEdgeNodes(VertNodes, ElementFaces, FaceNodes);
702  for (int j = 0; j < Nodes.num_elements(); ++j)
703  {
704  nodeList.push_back(VertNodes[Nodes[j]]);
705  }
706 
707  // Create element
708  ElmtConfig conf(elType, 1, true, true);
709  ElementSharedPtr E =
710  GetElementFactory().CreateInstance(elType, conf, nodeList, tags);
711 
712  m_mesh->m_element[E->GetDim()].push_back(E);
713 }
714 
715 void InputTec::GenElement3D(vector<NodeSharedPtr> &VertNodes,
716  int i,
717  vector<int> &ElementFaces,
718  vector<vector<int> > &FaceNodes,
719  int nComposite,
720  bool DoOrient)
721 {
723  // set up Node list
724  Array<OneD, int> Nodes = SortFaceNodes(VertNodes, ElementFaces, FaceNodes);
725  int nnodes = Nodes.num_elements();
726  map<LibUtilities::ShapeType, int> domainComposite;
727 
728  // Set Nodes -- Not sure we need this so could
729  // m_mesh->m_node = VertNodes;
730 
731  // element type
732  if (nnodes == 4)
733  {
735  }
736  else if (nnodes == 5)
737  {
738  elType = LibUtilities::ePyramid;
739  }
740  else if (nnodes == 6)
741  {
742  elType = LibUtilities::ePrism;
743  }
744  else
745  {
746 
747  ASSERTL0(false, "Not set up for elements which are not Tets or Prism");
748  }
749 
750  // Create element tags
751  vector<int> tags;
752  tags.push_back(nComposite);
753 
754  // make unique node list
755  vector<NodeSharedPtr> nodeList;
756  for (int j = 0; j < Nodes.num_elements(); ++j)
757  {
758  nodeList.push_back(VertNodes[Nodes[j]]);
759  }
760 
761  // Create element
762  if (elType != LibUtilities::ePyramid)
763  {
764  ElmtConfig conf(elType, 1, true, true, DoOrient);
765  ElementSharedPtr E =
766  GetElementFactory().CreateInstance(elType, conf, nodeList, tags);
767 
768  m_mesh->m_element[E->GetDim()].push_back(E);
769  }
770  else
771  {
772  cout << "Warning: Pyramid detected " << endl;
773  }
774 }
775 
776 Array<OneD, int> InputTec::SortEdgeNodes(vector<NodeSharedPtr> &Vnodes,
777  vector<int> &ElementFaces,
778  vector<vector<int> > &FaceNodes)
779 {
780  int i, j;
781  Array<OneD, int> returnval;
782 
783  if (ElementFaces.size() == 3) // Triangle
784  {
785  returnval = Array<OneD, int>(3);
786 
787  returnval[0] = FaceNodes[ElementFaces[0]][0];
788  returnval[1] = FaceNodes[ElementFaces[0]][1];
789 
790  // Find third node index;
791  for (i = 0; i < 2; ++i)
792  {
793  if ((FaceNodes[ElementFaces[1]][i] != returnval[0]) &&
794  (FaceNodes[ElementFaces[1]][i] != returnval[1]))
795  {
796  returnval[2] = FaceNodes[ElementFaces[1]][i];
797  break;
798  }
799  }
800  }
801  else if (ElementFaces.size() == 4) // quadrilateral
802  {
803  returnval = Array<OneD, int>(4);
804 
805  int indx0 = FaceNodes[ElementFaces[0]][0];
806  int indx1 = FaceNodes[ElementFaces[0]][1];
807  int indx2, indx3;
808 
809  indx2 = indx3 = -1;
810  // Find third, fourth node index;
811  for (j = 1; j < 4; ++j)
812  {
813  for (i = 0; i < 2; ++i)
814  {
815  if ((FaceNodes[ElementFaces[j]][i] != indx0) &&
816  (FaceNodes[ElementFaces[j]][i] != indx1))
817  {
818  if (indx2 == -1)
819  {
820  indx2 = FaceNodes[ElementFaces[j]][i];
821  }
822  else if (indx2 != -1)
823  {
824  if (FaceNodes[ElementFaces[j]][i] != indx2)
825  {
826  indx3 = FaceNodes[ElementFaces[j]][i];
827  }
828  }
829  }
830  }
831  }
832 
833  ASSERTL1((indx2 != -1) && (indx3 != -1),
834  "Failed to find vertex 3 or 4");
835 
836  // calculate 0-1,
837  Node a = *(Vnodes[indx1]) - *(Vnodes[indx0]);
838  // calculate 0-2,
839  Node b = *(Vnodes[indx2]) - *(Vnodes[indx0]);
840  Node acurlb = a.curl(b);
841 
842  // calculate 2-1,
843  Node c = *(Vnodes[indx1]) - *(Vnodes[indx2]);
844  // calculate 3-2,
845  Node d = *(Vnodes[indx3]) - *(Vnodes[indx2]);
846  Node acurld = a.curl(d);
847 
848  NekDouble acurlb_dot_acurld = acurlb.dot(acurld);
849  if (acurlb_dot_acurld > 0.0)
850  {
851  returnval[0] = indx0;
852  returnval[1] = indx1;
853  returnval[2] = indx2;
854  returnval[3] = indx3;
855  }
856  else
857  {
858  returnval[0] = indx0;
859  returnval[1] = indx1;
860  returnval[2] = indx3;
861  returnval[3] = indx2;
862  }
863  }
864 
865  return returnval;
866 }
867 
868 Array<OneD, int> InputTec::SortFaceNodes(vector<NodeSharedPtr> &Vnodes,
869  vector<int> &ElementFaces,
870  vector<vector<int> > &FaceNodes)
871 {
872 
873  int i, j;
874  Array<OneD, int> returnval;
875 
876  if (ElementFaces.size() == 4) // Tetrahedron
877  {
878  ASSERTL1(FaceNodes[ElementFaces[0]].size() == 3,
879  "Face is not triangular");
880 
881  returnval = Array<OneD, int>(4);
882 
883  int indx0 = FaceNodes[ElementFaces[0]][0];
884  int indx1 = FaceNodes[ElementFaces[0]][1];
885  int indx2 = FaceNodes[ElementFaces[0]][2];
886  int indx3 = -1;
887 
888  // calculate 0-1,
889  Node a = *(Vnodes[indx1]) - *(Vnodes[indx0]);
890  // calculate 0-2,
891  Node b = *(Vnodes[indx2]) - *(Vnodes[indx0]);
892 
893  // Find fourth node index;
894  ASSERTL1(FaceNodes[ElementFaces[1]].size() == 3,
895  "Face is not triangular");
896  for (i = 0; i < 3; ++i)
897  {
898 
899  if ((FaceNodes[ElementFaces[1]][i] != indx0) &&
900  (FaceNodes[ElementFaces[1]][i] != indx1) &&
901  (FaceNodes[ElementFaces[1]][i] != indx2))
902  {
903  indx3 = FaceNodes[ElementFaces[1]][i];
904  break;
905  }
906  }
907 
908  // calculate 0-3,
909  Node c = *(Vnodes[indx3]) - *(Vnodes[indx0]);
910  Node acurlb = a.curl(b);
911 
912  NekDouble acurlb_dotc = acurlb.dot(c);
913  if (acurlb_dotc < 0.0)
914  {
915  returnval[0] = indx0;
916  returnval[1] = indx1;
917  returnval[2] = indx2;
918  returnval[3] = indx3;
919  }
920  else
921  {
922  returnval[0] = indx1;
923  returnval[1] = indx0;
924  returnval[2] = indx2;
925  returnval[3] = indx3;
926  }
927  }
928  else if (ElementFaces.size() == 5) // prism or pyramid
929  {
930  int triface0, triface1;
931  int quadface0, quadface1, quadface2;
932  bool isPrism = true;
933 
934  // find ids of tri faces and first quad face
935  triface0 = triface1 = -1;
936  quadface0 = quadface1 = quadface2 = -1;
937  for (i = 0; i < 5; ++i)
938  {
939  if (FaceNodes[ElementFaces[i]].size() == 3)
940  {
941  if (triface0 == -1)
942  {
943  triface0 = i;
944  }
945  else if (triface1 == -1)
946  {
947  triface1 = i;
948  }
949  else
950  {
951  isPrism = false;
952  }
953  }
954 
955  if (FaceNodes[ElementFaces[i]].size() == 4)
956  {
957  if (quadface0 == -1)
958  {
959  quadface0 = i;
960  }
961  else if (quadface1 == -1)
962  {
963  quadface1 = i;
964  }
965  else if (quadface2 == -1)
966  {
967  quadface2 = i;
968  }
969  }
970  }
971 
972  if (isPrism) // Prism
973  {
974  returnval = Array<OneD, int>(6);
975  }
976  else // Pyramid
977  {
978  returnval = Array<OneD, int>(5);
979  }
980 
981  // find matching nodes between triface0 and triquad0
982  int indx0, indx1, indx2, indx3, indx4;
983 
984  indx0 = indx1 = indx2 = indx3 = indx4 = -1;
985  // Loop over all quad nodes and if they match any
986  // triangular nodes If they do set these to indx0 and
987  // indx1 and if not set it to indx2, indx3
988 
989  for (i = 0; i < 4; ++i)
990  {
991  for (j = 0; j < 3; ++j)
992  {
993  if (FaceNodes[ElementFaces[triface0]][j] ==
994  FaceNodes[ElementFaces[quadface0]][i])
995  {
996  break; // same node break
997  }
998  }
999 
1000  if (j == 3) // Vertex not in quad face
1001  {
1002  if (indx2 == -1)
1003  {
1004  indx2 = FaceNodes[ElementFaces[quadface0]][i];
1005  }
1006  else if (indx3 == -1)
1007  {
1008  indx3 = FaceNodes[ElementFaces[quadface0]][i];
1009  }
1010  else
1011  {
1012  ASSERTL0(
1013  false,
1014  "More than two vertices do not match triangular face");
1015  }
1016  }
1017  else // if found match then set indx0,indx1;
1018  {
1019  if (indx0 == -1)
1020  {
1021  indx0 = FaceNodes[ElementFaces[quadface0]][i];
1022  }
1023  else
1024  {
1025  indx1 = FaceNodes[ElementFaces[quadface0]][i];
1026  }
1027  }
1028  }
1029 
1030  // Finally check for top vertex
1031  for (int i = 0; i < 3; ++i)
1032  {
1033  if ((FaceNodes[ElementFaces[triface0]][i] != indx0) &&
1034  (FaceNodes[ElementFaces[triface0]][i] != indx1) &&
1035  (FaceNodes[ElementFaces[triface0]][i] != indx2))
1036  {
1037  indx4 = FaceNodes[ElementFaces[triface0]][i];
1038  break;
1039  }
1040  }
1041 
1042  // calculate 0-1,
1043  Node a = *(Vnodes[indx1]) - *(Vnodes[indx0]);
1044  // calculate 0-4,
1045  Node b = *(Vnodes[indx4]) - *(Vnodes[indx0]);
1046  // calculate 0-2,
1047  Node c = *(Vnodes[indx2]) - *(Vnodes[indx0]);
1048  Node acurlb = a.curl(b);
1049 
1050  NekDouble acurlb_dotc = acurlb.dot(c);
1051  if (acurlb_dotc < 0.0)
1052  {
1053  returnval[0] = indx0;
1054  returnval[1] = indx1;
1055  returnval[4] = indx4;
1056  }
1057  else
1058  {
1059  returnval[0] = indx1;
1060  returnval[1] = indx0;
1061  returnval[4] = indx4;
1062  }
1063 
1064  // check to see if two vertices are shared between one of the other
1065  // faces
1066  // to define which is indx2 and indx3
1067 
1068  int cnt = 0;
1069  for (int i = 0; i < 4; ++i)
1070  {
1071  if ((FaceNodes[ElementFaces[quadface1]][i] == returnval[1]) ||
1072  (FaceNodes[ElementFaces[quadface1]][i] == indx2))
1073  {
1074  cnt++;
1075  }
1076  }
1077 
1078  if (cnt == 2) // have two matching vertices
1079  {
1080  returnval[2] = indx2;
1081  returnval[3] = indx3;
1082  }
1083  else
1084  {
1085  cnt = 0;
1086  for (int i = 0; i < 4; ++i)
1087  {
1088  if ((FaceNodes[ElementFaces[quadface2]][i] == returnval[1]) ||
1089  (FaceNodes[ElementFaces[quadface2]][i] == indx2))
1090  {
1091  cnt++;
1092  }
1093  }
1094 
1095  if (cnt != 2) // neither of the other faces has two matching nodes
1096  // so reverse
1097  {
1098  returnval[2] = indx3;
1099  returnval[3] = indx2;
1100  }
1101  else // have two matching vertices
1102  {
1103  returnval[2] = indx2;
1104  returnval[3] = indx3;
1105  }
1106  }
1107 
1108  if (isPrism == true)
1109  {
1110  // finally need to find last vertex from second triangular face.
1111  for (int i = 0; i < 3; ++i)
1112  {
1113  if ((FaceNodes[ElementFaces[triface1]][i] != indx2) &&
1114  (FaceNodes[ElementFaces[triface1]][i] != indx3) &&
1115  (FaceNodes[ElementFaces[triface1]][i] != indx3))
1116  {
1117  returnval[5] = FaceNodes[ElementFaces[triface1]][i];
1118  break;
1119  }
1120  }
1121  }
1122  }
1123  else
1124  {
1125  ASSERTL0(false, "SortFaceNodes not set up for this number of faces");
1126  }
1127 
1128  return returnval;
1129 }
1130 }
1131 }
NEKMESHUTILS_EXPORT NekDouble dot(const Node &pSrc) const
Definition: Node.h:153
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
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
void ResetNodes(vector< NodeSharedPtr > &Nodes, Array< OneD, vector< int > > &ElementFaces, vector< vector< int > > &FaceNodes)
STL namespace.
MeshSharedPtr m_mesh
Mesh object.
Array< OneD, int > SortFaceNodes(vector< NodeSharedPtr > &Nodes, vector< int > &ElementFaces, vector< vector< int > > &FaceNodes)
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.
virtual void Process()
Populate and validate required data structures.
virtual void ProcessVertices()
Extract element vertices.
virtual void ProcessElements()
Generate element IDs.
void GenElement3D(vector< NodeSharedPtr > &Nodes, int i, vector< int > &ElementFaces, vector< vector< int > > &FaceNodes, int ncomposite, bool DoOrient)
NEKMESHUTILS_EXPORT Node curl(const Node &pSrc) const
Definition: Node.h:158
double NekDouble
Array< OneD, int > SortEdgeNodes(vector< NodeSharedPtr > &Nodes, vector< int > &ElementFaces, vector< vector< int > > &FaceNodes)
virtual void ProcessComposites()
Generate composites.
void PrintSummary()
Print summary of elements.
void ReadZone(int &nComposite)
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
static void PrismLineFaces(int prismid, map< int, int > &facelist, vector< vector< int > > &FacesToPrisms, vector< vector< int > > &PrismsToFaces, vector< bool > &PrismDone)
Definition: InputStar.cpp:422
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Definition: Mesh.h:137
virtual void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:52
void GenElement2D(vector< NodeSharedPtr > &Nodes, int i, vector< int > &ElementFaces, vector< vector< int > > &FaceNodes, int ncomposite)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
void OpenStream()
Open a file for input.
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