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