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