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