Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 << "InputCCM: 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, triface2, triface3;
676  int quadface0, quadface1, quadface2;
677  bool isPrism = true;
678 
679  // find ids of tri faces and first quad face
680  triface0 = triface1 = triface2 = triface3 = -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 if (triface2 == -1)
695  {
696  triface2 = i;
697  isPrism = false;
698  }
699  else if (triface3 == -1)
700  {
701  triface3 = i;
702  }
703  }
704 
705  if (FaceNodes[ElementFaces[i]].size() == 4)
706  {
707  if (quadface0 == -1)
708  {
709  quadface0 = i;
710  }
711  else if (quadface1 == -1)
712  {
713  quadface1 = i;
714  }
715  else if (quadface2 == -1)
716  {
717  quadface2 = i;
718  }
719  }
720  }
721 
722  if (isPrism) // Prism
723  {
724  returnval = Array<OneD, int>(6);
725  ASSERTL1(quadface0 != -1, "Quad face 0 not found");
726  ASSERTL1(quadface1 != -1, "Quad face 1 not found");
727  ASSERTL1(quadface2 != -1, "Quad face 2 not found");
728  ASSERTL1(triface0 != -1, "Tri face 0 not found");
729  ASSERTL1(triface1 != -1, "Tri face 1 not found");
730  }
731  else // Pyramid
732  {
733  ASSERTL1(quadface0 != -1, "Quad face 0 not found");
734  ASSERTL1(triface0 != -1, "Tri face 0 not found");
735  ASSERTL1(triface1 != -1, "Tri face 1 not found");
736  ASSERTL1(triface2 != -1, "Tri face 2 not found");
737  ASSERTL1(triface3 != -1, "Tri face 3 not found");
738  ASSERTL0(false,"Pyramids still not sorted");
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  else if (indx3 == -1)
768  {
769  indx3 = FaceNodes[ElementFaces[quadface0]][i];
770  }
771  else
772  {
773  ASSERTL0(
774  false,
775  "More than two vertices do not match triangular face");
776  }
777  }
778  else // if found match then set indx0,indx1;
779  {
780  if (indx0 == -1)
781  {
782  indx0 = FaceNodes[ElementFaces[quadface0]][i];
783  }
784  else
785  {
786  indx1 = FaceNodes[ElementFaces[quadface0]][i];
787  }
788  }
789  }
790 
791  // Finally check for top vertex
792  for (int i = 0; i < 3; ++i)
793  {
794  if ((FaceNodes[ElementFaces[triface0]][i] != indx0) &&
795  (FaceNodes[ElementFaces[triface0]][i] != indx1) &&
796  (FaceNodes[ElementFaces[triface0]][i] != indx2))
797  {
798  indx4 = FaceNodes[ElementFaces[triface0]][i];
799  break;
800  }
801  }
802 
803  // calculate 0-1,
804  Node a = *(Vnodes[indx1]) - *(Vnodes[indx0]);
805  // calculate 0-4,
806  Node b = *(Vnodes[indx4]) - *(Vnodes[indx0]);
807  // calculate 0-2,
808  Node c = *(Vnodes[indx2]) - *(Vnodes[indx0]);
809  Node acurlb = a.curl(b);
810 
811  NekDouble acurlb_dotc = acurlb.dot(c);
812  if (acurlb_dotc < 0.0)
813  {
814  returnval[0] = indx0;
815  returnval[1] = indx1;
816  returnval[4] = indx4;
817  }
818  else
819  {
820  returnval[0] = indx1;
821  returnval[1] = indx0;
822  returnval[4] = indx4;
823  }
824 
825  // check to see if two vertices are shared between one of the other
826  // faces
827  // to define which is indx2 and indx3
828 
829  int cnt = 0;
830  for (int i = 0; i < 4; ++i)
831  {
832  if ((FaceNodes[ElementFaces[quadface1]][i] == returnval[1]) ||
833  (FaceNodes[ElementFaces[quadface1]][i] == indx2))
834  {
835  cnt++;
836  }
837  }
838 
839  if (cnt == 2) // have two matching vertices
840  {
841  returnval[2] = indx2;
842  returnval[3] = indx3;
843  }
844  else
845  {
846  cnt = 0;
847  for (int i = 0; i < 4; ++i)
848  {
849  if ((FaceNodes[ElementFaces[quadface2]][i] == returnval[1]) ||
850  (FaceNodes[ElementFaces[quadface2]][i] == indx2))
851  {
852  cnt++;
853  }
854  }
855 
856  if (cnt != 2) // neither of the other faces has two matching nodes
857  // so reverse
858  {
859  returnval[2] = indx3;
860  returnval[3] = indx2;
861  }
862  else // have two matching vertices
863  {
864  returnval[2] = indx2;
865  returnval[3] = indx3;
866  }
867  }
868 
869  if (isPrism == true)
870  {
871  // finally need to find last vertex from second triangular face.
872  for (int i = 0; i < 3; ++i)
873  {
874  if ((FaceNodes[ElementFaces[triface1]][i] != indx2) &&
875  (FaceNodes[ElementFaces[triface1]][i] != indx3) &&
876  (FaceNodes[ElementFaces[triface1]][i] != indx3))
877  {
878  returnval[5] = FaceNodes[ElementFaces[triface1]][i];
879  break;
880  }
881  }
882  }
883  }
884  else
885  {
886  ASSERTL0(false, "SortFaceNodes not set up for this number of faces");
887  }
888 
889  return returnval;
890 }
891 
892 // initialise and read ccm file to ccm structure
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(
929  &m_ccmErr, m_ccmProcessor, &vertices, &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,
946  vertices,
947  &dims,
948  &scale,
949  &mapID,
950  verts,
951  CCMIOINDEXC(0),
952  CCMIOINDEXC(0 + nVertices));
953  CCMIOReadMap(
954  &m_ccmErr, mapID, mapData, CCMIOINDEXC(0), CCMIOINDEXC(0 + nVertices));
955 
956  for (int i = 0; i < nVertices; ++i)
957  {
958  Nodes.push_back(boost::shared_ptr<Node>(
959  new Node(i, verts[3 * i], verts[3 * i + 1], verts[3 * i + 2])));
960  }
961 }
962 
963 void InputStar::ReadInternalFaces(map<int, vector<int> > &FacesNodes,
964  Array<OneD, vector<int> > &ElementFaces)
965 {
966 
967  CCMIOID mapID, id;
968  CCMIOSize_t nFaces, size;
969  vector<int> faces, faceCells, mapData;
970 
971  // Read the internal faces.
972  CCMIOGetEntity(&m_ccmErr, m_ccmTopology, kCCMIOInternalFaces, 0, &id);
973  CCMIOEntitySize(&m_ccmErr, id, &nFaces, NULL);
974 
975  int nf = TOINT64(nFaces);
976  mapData.resize(nf);
977  faceCells.resize(2 * nf);
978 
979  CCMIOReadFaces(&m_ccmErr,
980  id,
981  kCCMIOInternalFaces,
982  NULL,
983  &size,
984  NULL,
985  CCMIOINDEXC(kCCMIOStart),
986  CCMIOINDEXC(kCCMIOEnd));
987  faces.resize(TOINT64(size));
988  CCMIOReadFaces(&m_ccmErr,
989  id,
990  kCCMIOInternalFaces,
991  &mapID,
992  NULL,
993  &faces[0],
994  CCMIOINDEXC(kCCMIOStart),
995  CCMIOINDEXC(kCCMIOEnd));
996  CCMIOReadFaceCells(&m_ccmErr,
997  id,
998  kCCMIOInternalFaces,
999  &faceCells[0],
1000  CCMIOINDEXC(kCCMIOStart),
1001  CCMIOINDEXC(kCCMIOEnd));
1002  CCMIOReadMap(&m_ccmErr,
1003  mapID,
1004  &mapData[0],
1005  CCMIOINDEXC(kCCMIOStart),
1006  CCMIOINDEXC(kCCMIOEnd));
1007 
1008  // Add face nodes
1009  int cnt = 0;
1010  for (int i = 0; i < nf; ++i)
1011  {
1012  vector<int> Fnodes;
1013  int j;
1014  if (cnt < faces.size())
1015  {
1016  int nv = faces[cnt];
1017  ASSERTL0(nv <= 4,
1018  "Can only handle meshes with "
1019  "up to four nodes per face");
1020 
1021  for (j = 0; j < nv; ++j)
1022  {
1023  if (cnt + 1 + j < faces.size())
1024  {
1025  Fnodes.push_back(faces[cnt + 1 + j] - 1);
1026  }
1027  }
1028  cnt += nv + 1;
1029  }
1030  FacesNodes[mapData[i] - 1] = Fnodes;
1031  }
1032 
1033  // find number of elements;
1034  int nelmt = 0;
1035  for (int i = 0; i < faceCells.size(); ++i)
1036  {
1037  nelmt = max(nelmt, faceCells[i]);
1038  }
1039 
1040  ElementFaces = Array<OneD, vector<int> >(nelmt);
1041  for (int i = 0; i < nf; ++i)
1042  {
1043  // left element
1044  if (faceCells[2 * i])
1045  {
1046  ElementFaces[faceCells[2 * i] - 1].push_back(mapData[i] - 1);
1047  }
1048 
1049  // right element
1050  if (faceCells[2 * i + 1])
1051  {
1052  ElementFaces[faceCells[2 * i + 1] - 1].push_back(mapData[i] - 1);
1053  }
1054  }
1055 }
1056 
1057 void InputStar::ReadBoundaryFaces(vector<vector<int> > &BndElementFaces,
1058  map<int, vector<int> > &FacesNodes,
1059  Array<OneD, vector<int> > &ElementFaces,
1060  vector<string> &Facelabels)
1061 {
1062  // Read the boundary faces.
1063  CCMIOSize_t index = CCMIOSIZEC(0);
1064  CCMIOID mapID, id;
1065  CCMIOSize_t nFaces, size;
1066  vector<int> faces, faceCells, mapData;
1067  vector<string> facelabel;
1068 
1069  while (CCMIONextEntity(
1070  NULL, m_ccmTopology, kCCMIOBoundaryFaces, &index, &id) ==
1071  kCCMIONoErr)
1072  {
1073  int boundaryVal;
1074 
1075  CCMIOEntitySize(&m_ccmErr, id, &nFaces, NULL);
1076  int nf = TOINT64(nFaces);
1077  mapData.resize(nf);
1078  faceCells.resize(nf);
1079  CCMIOReadFaces(&m_ccmErr,
1080  id,
1081  kCCMIOBoundaryFaces,
1082  NULL,
1083  &size,
1084  NULL,
1085  CCMIOINDEXC(kCCMIOStart),
1086  CCMIOINDEXC(kCCMIOEnd));
1087 
1088  faces.resize(TOINT64(size));
1089  CCMIOReadFaces(&m_ccmErr,
1090  id,
1091  kCCMIOBoundaryFaces,
1092  &mapID,
1093  NULL,
1094  &faces[0],
1095  CCMIOINDEXC(kCCMIOStart),
1096  CCMIOINDEXC(kCCMIOEnd));
1097  CCMIOReadFaceCells(&m_ccmErr,
1098  id,
1099  kCCMIOBoundaryFaces,
1100  &faceCells[0],
1101  CCMIOINDEXC(kCCMIOStart),
1102  CCMIOINDEXC(kCCMIOEnd));
1103  CCMIOReadMap(&m_ccmErr,
1104  mapID,
1105  &mapData[0],
1106  CCMIOINDEXC(kCCMIOStart),
1107  CCMIOINDEXC(kCCMIOEnd));
1108 
1109  CCMIOGetEntityIndex(&m_ccmErr, id, &boundaryVal);
1110 
1111  // check to see if we have a label for this boundary faces
1112  int size;
1113  char *name;
1114  if (CCMIOReadOptstr(NULL, id, "Label", &size, NULL) == kCCMIONoErr)
1115  {
1116  name = new char[size + 1];
1117  CCMIOReadOptstr(NULL, id, "Label", NULL, name);
1118  Facelabels.push_back(string(name));
1119  }
1120  else
1121  {
1122  Facelabels.push_back("Not known");
1123  }
1124 
1125  // Add face nodes
1126  int cnt = 0;
1127  for (int i = 0; i < nf; ++i)
1128  {
1129  vector<int> Fnodes;
1130  int j;
1131  if (cnt < faces.size())
1132  {
1133  int nv = faces[cnt];
1134  ASSERTL0(nv <= 4,
1135  "Can only handle meshes with "
1136  "up to four nodes per face");
1137 
1138  for (j = 0; j < nv; ++j)
1139  {
1140  if (cnt + 1 + j < faces.size())
1141  {
1142  Fnodes.push_back(faces[cnt + 1 + j] - 1);
1143  }
1144  }
1145  cnt += nv + 1;
1146  }
1147  FacesNodes[mapData[i] - 1] = Fnodes;
1148  }
1149 
1150  vector<int> BndFaces;
1151  for (int i = 0; i < nf; ++i)
1152  {
1153  if (faceCells[i])
1154  {
1155  ElementFaces[faceCells[i] - 1].push_back(mapData[i] - 1);
1156  }
1157  BndFaces.push_back(mapData[i] - 1);
1158  }
1159  BndElementFaces.push_back(BndFaces);
1160  }
1161 }
1162 }
1163 }
NEKMESHUTILS_EXPORT NekDouble dot(const Node &pSrc) const
Definition: Node.h:159
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Basic information about an element.
Definition: ElementConfig.h:50
virtual void Process()
Populate and validate required data structures.
Definition: InputStar.cpp:77
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
STL namespace.
pair< ModuleType, string > ModuleKey
void GenElement3D(std::vector< NekMeshUtils::NodeSharedPtr > &Nodes, int i, std::vector< int > &ElementFaces, std::map< int, std::vector< int > > &FaceNodes, int ncomposite, bool DoOrient)
Definition: InputStar.cpp:498
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
Represents a point in the domain.
Definition: Node.h:60
NEKMESHUTILS_EXPORT void PrintSummary()
Print summary of elements.
void ReadBoundaryFaces(std::vector< std::vector< int > > &BndElementFaces, std::map< int, std::vector< int > > &FacesNodes, Array< OneD, std::vector< int > > &ElementFaces, std::vector< std::string > &facelabels)
Definition: InputStar.cpp:1057
virtual NEKMESHUTILS_EXPORT void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
virtual NEKMESHUTILS_EXPORT void ProcessElements()
Generate element IDs.
Represents a command-line configuration option.
void ReadNodes(std::vector< NekMeshUtils::NodeSharedPtr > &Nodes)
Definition: InputStar.cpp:923
NEKMESHUTILS_EXPORT Node curl(const Node &pSrc) const
Definition: Node.h:164
double NekDouble
std::map< std::string, ConfigOption > m_config
List of configuration values.
void ReadInternalFaces(std::map< int, std::vector< int > > &FacesNodes, Array< OneD, std::vector< int > > &ElementFaces)
Definition: InputStar.cpp:963
string msg
print "Adding",units.name, hash(units), units.description(), print "(was",id(_u),"now",id(units),")" Ensure referenced units exist
Definition: pycml.py:3834
void GenElement2D(std::vector< NekMeshUtils::NodeSharedPtr > &Nodes, int i, std::vector< int > &FaceNodes, int ncomposite)
Definition: InputStar.cpp:458
Abstract base class for input modules.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
static void PrismLineFaces(int prismid, map< int, int > &facelist, vector< vector< int > > &FacesToPrisms, vector< vector< int > > &PrismsToFaces, vector< bool > &PrismDone)
Definition: InputStar.cpp:422
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Definition: Mesh.h:147
Array< OneD, int > SortFaceNodes(std::vector< NekMeshUtils::NodeSharedPtr > &Nodes, std::vector< int > &ElementFaces, std::map< int, std::vector< int > > &FaceNodes)
Definition: InputStar.cpp:613
virtual NEKMESHUTILS_EXPORT void ProcessVertices()
Extract element vertices.
Array< OneD, int > SortEdgeNodes(std::vector< NekMeshUtils::NodeSharedPtr > &Nodes, std::vector< int > &FaceNodes)
Definition: InputStar.cpp:559
virtual NEKMESHUTILS_EXPORT void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
static char const kDefaultState[]
Definition: InputStar.cpp:47
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
void ResetNodes(std::vector< NekMeshUtils::NodeSharedPtr > &Nodes, Array< OneD, std::vector< int > > &ElementFaces, std::map< int, std::vector< int > > &FaceNodes)
Definition: InputStar.cpp:210
std::pair< ModuleType, std::string > ModuleKey
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
virtual NEKMESHUTILS_EXPORT void ProcessComposites()
Generate composites.
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215