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