36 #include <boost/algorithm/string.hpp>
54 "Reads Tecplot polyhedron ascii format converted from Star CCM (.dat).");
81 cout <<
"InputStarTec: Start reading file..." << endl;
95 if (line.find(
"ZONE") != string::npos)
105 if (line.find(
"ZONE") != string::npos)
124 int nfaces, nnodes, nelements;
131 nnodes = nfaces = nelements = 0;
136 boost::to_upper(line);
139 if (sscanf(line.c_str(),
"%lf", &value) == 1)
144 if ((line.find(
"NODES") != string::npos) &&
145 (line.find(
"TOTALNUMFACENODES") == string::npos))
151 start = tag.find(
"NODES=");
152 end = tag.find_first_of(
',', start);
153 nnodes = atoi(tag.substr(start + 6, end).c_str());
156 if ((line.find(
"FACES") != string::npos) &&
157 (line.find(
"NUMCONNECTEDBOUNDARYFACES") == string::npos))
163 start = tag.find(
"FACES=");
164 end = tag.find_first_of(
',', start);
165 nfaces = atoi(tag.substr(start + 6, end).c_str());
168 if (line.find(
"ELEMENTS") != string::npos)
174 start = tag.find(
"ELEMENTS=");
175 end = tag.find_first_of(
',', start);
176 nelements = atoi(tag.substr(start + 9, end).c_str());
179 if (line.find(
"ZONETYPE") != string::npos)
184 if ((line.find(
"FEPOLYGON") == string::npos) &&
185 (line.find(
"FEPOLYHEDRON") == string::npos))
188 "Routine only set up for FEPolygon or FEPolyhedron");
197 cout <<
"Setting up zone " << zcnt++;
199 int nodeCount = 3 * nnodes;
200 vector<NekDouble> nodeLocs;
202 while (nodeCount > 0 && !
m_mshFile.eof())
208 nodeLocs.push_back(value);
217 ASSERTL0(nodeLocs.size() == 3*nnodes,
"Unable to read correct number of "
218 "nodes from Tecplot file");
220 std::vector<NodeSharedPtr> Nodes;
221 for (i = 0; i < nnodes; ++i)
224 boost::shared_ptr<Node>(
225 new Node(i, nodeLocs[i], nodeLocs[i+nnodes],
226 nodeLocs[i+2*nnodes])));
231 if (line.find(
"node count per face") == string::npos)
233 if (line.find(
"face nodes") == string::npos)
242 vector<int> Nodes_per_face;
243 if (line.find(
"node count per face") != string::npos)
246 for (i = 0; i < nfaces; ++i)
250 "Can only handle meshes with "
251 "up to four nodes per face");
252 Nodes_per_face.push_back(nodes);
259 if (line.find(
"face nodes") == string::npos)
266 vector<vector<int> > FaceNodes;
268 if (line.find(
"face nodes") != string::npos)
271 for (i = 0; i < nfaces; ++i)
275 int nodes = (Nodes_per_face.size()) ? Nodes_per_face[i] : 2;
279 for (
int j = 0; j < nodes; ++j)
284 Fnodes.push_back(nodeID - 1);
287 FaceNodes.push_back(Fnodes);
292 ASSERTL0(
false,
"Failed to find face node section");
300 if (line.find(
"left elements") == string::npos)
305 if (line.find(
"left elements") != string::npos)
309 for (i = 0; i < nfaces; ++i)
315 ElementFaces[elmtID - 1].push_back(i);
321 ASSERTL0(
false,
"Left element not found");
326 if (line.find(
"right elements") == string::npos)
331 if (line.find(
"right elements") != string::npos)
336 for (i = 0; i < nfaces; ++i)
342 ElementFaces[elmtID - 1].push_back(i);
351 ASSERTL0(
false,
"Left element not found");
354 if (Nodes_per_face.size())
356 cout <<
" (3D) " << endl;
365 for (i = 0; i < nelements; ++i)
367 if (ElementFaces[i].size() > 4)
370 Nodes, i, ElementFaces[i], FaceNodes, nComposite,
true);
377 for (i = 0; i < nelements; ++i)
379 if (ElementFaces[i].size() == 4)
382 Nodes, i, ElementFaces[i], FaceNodes, nComposite,
true);
391 cout <<
" (2D)" << endl;
395 for (i = 0; i < Nodes.size(); ++i)
399 if (it ==
m_mesh->m_vertexSet.end())
401 ASSERTL0(
false,
"Failed to find face vertex in 3D list");
409 for (i = 0; i < nelements; ++i)
411 GenElement2D(Nodes, i, ElementFaces[i], FaceNodes, nComposite);
419 map<int, int> &facelist,
420 vector<vector<int> > &FacesToPrisms,
421 vector<vector<int> > &PrismsToFaces,
422 vector<bool> &PrismDone);
426 vector<vector<int> > &FaceNodes)
430 int face1_map[3] = {0, 1, 4};
431 int face3_map[3] = {3, 2, 5};
433 map<int, bool> FacesRenumbered;
436 vector<vector<int> > FaceToPrisms(FaceNodes.size());
437 vector<vector<int> > PrismToFaces(ElementFaces.num_elements());
438 map<int, int> Prisms;
443 for (i = 0; i < ElementFaces.num_elements(); ++i)
446 if (ElementFaces[i].size() == 5)
448 vector<int> LocTriFaces;
450 for (j = 0; j < ElementFaces[i].size(); ++j)
452 if (FaceNodes[ElementFaces[i][j]].size() == 3)
454 LocTriFaces.push_back(j);
458 if (LocTriFaces.size() == 2)
462 PrismToFaces[i].push_back(ElementFaces[i][LocTriFaces[0]]);
463 PrismToFaces[i].push_back(ElementFaces[i][LocTriFaces[1]]);
465 FaceToPrisms[ElementFaces[i][LocTriFaces[0]]].push_back(i);
466 FaceToPrisms[ElementFaces[i][LocTriFaces[1]]].push_back(i);
471 vector<bool> FacesDone(FaceNodes.size(),
false);
472 vector<bool> PrismDone(ElementFaces.num_elements(),
false);
477 for (PrismIt = Prisms.begin(); PrismIt != Prisms.end(); ++PrismIt)
479 int elmtid = PrismIt->first;
480 map<int, int> facelist;
483 if (PrismDone[elmtid])
491 elmtid, facelist, FaceToPrisms, PrismToFaces, PrismDone);
494 for (faceIt = facelist.begin(); faceIt != facelist.end(); faceIt++)
496 int faceid = faceIt->second;
498 for (i = 0; i < FaceToPrisms[faceid].size(); ++i)
500 int prismid = FaceToPrisms[faceid][i];
502 if ((FacesDone[PrismToFaces[prismid][0]] ==
true) &&
503 (FacesDone[PrismToFaces[prismid][1]] ==
true))
511 if ((FacesDone[PrismToFaces[prismid][0]] ==
false) &&
512 (FacesDone[PrismToFaces[prismid][1]] ==
false))
516 for (i = 0; i < 3; ++i)
518 if (NodeReordering[Nodes[face1_map[i]]] == -1)
520 NodeReordering[Nodes[face1_map[i]]] = nodeid++;
524 for (i = 0; i < 3; ++i)
526 if (NodeReordering[Nodes[face3_map[i]]] == -1)
528 NodeReordering[Nodes[face3_map[i]]] = nodeid++;
532 else if ((FacesDone[PrismToFaces[prismid][0]] ==
false) &&
533 (FacesDone[PrismToFaces[prismid][1]] ==
true))
536 int max_id1, max_id2;
538 max_id1 = (NodeReordering[Nodes[face3_map[0]]] <
539 NodeReordering[Nodes[face3_map[1]]])
542 max_id2 = (NodeReordering[Nodes[face3_map[max_id1]]] <
543 NodeReordering[Nodes[face3_map[2]]])
548 int id0 = (max_id1 == 1) ? 0 : 1;
550 if (NodeReordering[Nodes[face1_map[id0]]] == -1)
552 NodeReordering[Nodes[face1_map[id0]]] = nodeid++;
555 if (NodeReordering[Nodes[face1_map[max_id1]]] == -1)
557 NodeReordering[Nodes[face1_map[max_id1]]] =
561 if (NodeReordering[Nodes[face1_map[max_id2]]] == -1)
563 NodeReordering[Nodes[face1_map[max_id2]]] =
567 else if ((FacesDone[PrismToFaces[prismid][0]] ==
true) &&
568 (FacesDone[PrismToFaces[prismid][1]] ==
false))
571 int max_id1, max_id2;
573 max_id1 = (NodeReordering[Nodes[face1_map[0]]] <
574 NodeReordering[Nodes[face1_map[1]]])
577 max_id2 = (NodeReordering[Nodes[face1_map[max_id1]]] <
578 NodeReordering[Nodes[face1_map[2]]])
583 int id0 = (max_id1 == 1) ? 0 : 1;
585 if (NodeReordering[Nodes[face3_map[id0]]] == -1)
587 NodeReordering[Nodes[face3_map[id0]]] = nodeid++;
590 if (NodeReordering[Nodes[face3_map[max_id1]]] == -1)
592 NodeReordering[Nodes[face3_map[max_id1]]] =
596 if (NodeReordering[Nodes[face3_map[max_id2]]] == -1)
598 NodeReordering[Nodes[face3_map[max_id2]]] =
608 for (i = 0; i < NodeReordering.num_elements(); ++i)
610 if (NodeReordering[i] == -1)
612 NodeReordering[i] = nodeid++;
616 ASSERTL1(nodeid == NodeReordering.num_elements(),
617 "Have not renumbered all nodes");
620 for (i = 0; i < FaceNodes.size(); ++i)
622 for (j = 0; j < FaceNodes[i].size(); ++j)
624 FaceNodes[i][j] = NodeReordering[FaceNodes[i][j]];
628 vector<NodeSharedPtr> save(Vnodes);
629 for (i = 0; i < Vnodes.size(); ++i)
631 Vnodes[NodeReordering[i]] = save[i];
632 Vnodes[NodeReordering[i]]->SetID(NodeReordering[i]);
637 map<int, int> &facelist,
638 vector<vector<int> > &FaceToPrisms,
639 vector<vector<int> > &PrismToFaces,
640 vector<bool> &PrismDone)
642 if (PrismDone[prismid] ==
false)
644 PrismDone[prismid] =
true;
647 int face = PrismToFaces[prismid][0];
648 facelist[face] = face;
649 for (
int i = 0; i < FaceToPrisms[face].size(); ++i)
659 face = PrismToFaces[prismid][1];
660 facelist[face] = face;
661 for (
int i = 0; i < FaceToPrisms[face].size(); ++i)
674 vector<int> &ElementFaces,
675 vector<vector<int> > &FaceNodes,
681 if (ElementFaces.size() == 3)
685 else if (ElementFaces.size() == 4)
691 ASSERTL0(
false,
"Not set up for elements which are not Tets or Prism");
696 tags.push_back(nComposite);
699 vector<NodeSharedPtr> nodeList;
701 for (
int j = 0; j < Nodes.num_elements(); ++j)
703 nodeList.push_back(VertNodes[Nodes[j]]);
711 m_mesh->m_element[E->GetDim()].push_back(E);
716 vector<int> &ElementFaces,
717 vector<vector<int> > &FaceNodes,
724 int nnodes = Nodes.num_elements();
725 map<LibUtilities::ShapeType, int> domainComposite;
735 else if (nnodes == 5)
739 else if (nnodes == 6)
746 ASSERTL0(
false,
"Not set up for elements which are not Tets or Prism");
751 tags.push_back(nComposite);
754 vector<NodeSharedPtr> nodeList;
755 for (
int j = 0; j < Nodes.num_elements(); ++j)
757 nodeList.push_back(VertNodes[Nodes[j]]);
763 ElmtConfig conf(elType, 1,
true,
true, DoOrient);
767 m_mesh->m_element[E->GetDim()].push_back(E);
771 cout <<
"Warning: Pyramid detected " << endl;
776 vector<int> &ElementFaces,
777 vector<vector<int> > &FaceNodes)
782 if (ElementFaces.size() == 3)
786 returnval[0] = FaceNodes[ElementFaces[0]][0];
787 returnval[1] = FaceNodes[ElementFaces[0]][1];
790 for (i = 0; i < 2; ++i)
792 if ((FaceNodes[ElementFaces[1]][i] != returnval[0]) &&
793 (FaceNodes[ElementFaces[1]][i] != returnval[1]))
795 returnval[2] = FaceNodes[ElementFaces[1]][i];
800 else if (ElementFaces.size() == 4)
804 int indx0 = FaceNodes[ElementFaces[0]][0];
805 int indx1 = FaceNodes[ElementFaces[0]][1];
810 for (j = 1; j < 4; ++j)
812 for (i = 0; i < 2; ++i)
814 if ((FaceNodes[ElementFaces[j]][i] != indx0) &&
815 (FaceNodes[ElementFaces[j]][i] != indx1))
819 indx2 = FaceNodes[ElementFaces[j]][i];
821 else if (indx2 != -1)
823 if (FaceNodes[ElementFaces[j]][i] != indx2)
825 indx3 = FaceNodes[ElementFaces[j]][i];
832 ASSERTL1((indx2 != -1) && (indx3 != -1),
833 "Failed to find vertex 3 or 4");
836 Node a = *(Vnodes[indx1]) - *(Vnodes[indx0]);
838 Node b = *(Vnodes[indx2]) - *(Vnodes[indx0]);
842 Node c = *(Vnodes[indx1]) - *(Vnodes[indx2]);
844 Node d = *(Vnodes[indx3]) - *(Vnodes[indx2]);
848 if (acurlb_dot_acurld > 0.0)
850 returnval[0] = indx0;
851 returnval[1] = indx1;
852 returnval[2] = indx2;
853 returnval[3] = indx3;
857 returnval[0] = indx0;
858 returnval[1] = indx1;
859 returnval[2] = indx3;
860 returnval[3] = indx2;
868 vector<int> &ElementFaces,
869 vector<vector<int> > &FaceNodes)
875 if (ElementFaces.size() == 4)
877 ASSERTL1(FaceNodes[ElementFaces[0]].size() == 3,
878 "Face is not triangular");
882 int indx0 = FaceNodes[ElementFaces[0]][0];
883 int indx1 = FaceNodes[ElementFaces[0]][1];
884 int indx2 = FaceNodes[ElementFaces[0]][2];
888 Node a = *(Vnodes[indx1]) - *(Vnodes[indx0]);
890 Node b = *(Vnodes[indx2]) - *(Vnodes[indx0]);
893 ASSERTL1(FaceNodes[ElementFaces[1]].size() == 3,
894 "Face is not triangular");
895 for (i = 0; i < 3; ++i)
898 if ((FaceNodes[ElementFaces[1]][i] != indx0) &&
899 (FaceNodes[ElementFaces[1]][i] != indx1) &&
900 (FaceNodes[ElementFaces[1]][i] != indx2))
902 indx3 = FaceNodes[ElementFaces[1]][i];
908 Node c = *(Vnodes[indx3]) - *(Vnodes[indx0]);
912 if (acurlb_dotc < 0.0)
914 returnval[0] = indx0;
915 returnval[1] = indx1;
916 returnval[2] = indx2;
917 returnval[3] = indx3;
921 returnval[0] = indx1;
922 returnval[1] = indx0;
923 returnval[2] = indx2;
924 returnval[3] = indx3;
927 else if (ElementFaces.size() == 5)
929 int triface0, triface1;
930 int quadface0, quadface1, quadface2;
934 triface0 = triface1 = -1;
935 quadface0 = quadface1 = quadface2 = -1;
936 for (i = 0; i < 5; ++i)
938 if (FaceNodes[ElementFaces[i]].size() == 3)
944 else if (triface1 == -1)
954 if (FaceNodes[ElementFaces[i]].size() == 4)
960 else if (quadface1 == -1)
964 else if (quadface2 == -1)
981 int indx0, indx1, indx2, indx3, indx4;
983 indx0 = indx1 = indx2 = indx3 = indx4 = -1;
988 for (i = 0; i < 4; ++i)
990 for (j = 0; j < 3; ++j)
992 if (FaceNodes[ElementFaces[triface0]][j] ==
993 FaceNodes[ElementFaces[quadface0]][i])
1003 indx2 = FaceNodes[ElementFaces[quadface0]][i];
1005 else if (indx3 == -1)
1007 indx3 = FaceNodes[ElementFaces[quadface0]][i];
1013 "More than two vertices do not match triangular face");
1020 indx0 = FaceNodes[ElementFaces[quadface0]][i];
1024 indx1 = FaceNodes[ElementFaces[quadface0]][i];
1030 for (
int i = 0; i < 3; ++i)
1032 if ((FaceNodes[ElementFaces[triface0]][i] != indx0) &&
1033 (FaceNodes[ElementFaces[triface0]][i] != indx1) &&
1034 (FaceNodes[ElementFaces[triface0]][i] != indx2))
1036 indx4 = FaceNodes[ElementFaces[triface0]][i];
1042 Node a = *(Vnodes[indx1]) - *(Vnodes[indx0]);
1044 Node b = *(Vnodes[indx4]) - *(Vnodes[indx0]);
1046 Node c = *(Vnodes[indx2]) - *(Vnodes[indx0]);
1050 if (acurlb_dotc < 0.0)
1052 returnval[0] = indx0;
1053 returnval[1] = indx1;
1054 returnval[4] = indx4;
1058 returnval[0] = indx1;
1059 returnval[1] = indx0;
1060 returnval[4] = indx4;
1068 for (
int i = 0; i < 4; ++i)
1070 if ((FaceNodes[ElementFaces[quadface1]][i] == returnval[1]) ||
1071 (FaceNodes[ElementFaces[quadface1]][i] == indx2))
1079 returnval[2] = indx2;
1080 returnval[3] = indx3;
1085 for (
int i = 0; i < 4; ++i)
1087 if ((FaceNodes[ElementFaces[quadface2]][i] == returnval[1]) ||
1088 (FaceNodes[ElementFaces[quadface2]][i] == indx2))
1097 returnval[2] = indx3;
1098 returnval[3] = indx2;
1102 returnval[2] = indx2;
1103 returnval[3] = indx3;
1107 if (isPrism ==
true)
1110 for (
int i = 0; i < 3; ++i)
1112 if ((FaceNodes[ElementFaces[triface1]][i] != indx2) &&
1113 (FaceNodes[ElementFaces[triface1]][i] != indx3) &&
1114 (FaceNodes[ElementFaces[triface1]][i] != indx3))
1116 returnval[5] = FaceNodes[ElementFaces[triface1]][i];
1124 ASSERTL0(
false,
"SortFaceNodes not set up for this number of faces");
NEKMESHUTILS_EXPORT NekDouble dot(const Node &pSrc) const
#define ASSERTL0(condition, msg)
Basic information about an element.
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.
pair< ModuleType, string > ModuleKey
MeshSharedPtr m_mesh
Mesh object.
ElementFactory & GetElementFactory()
Represents a point in the domain.
virtual NEKMESHUTILS_EXPORT void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
virtual NEKMESHUTILS_EXPORT void ProcessElements()
Generate element IDs.
NEKMESHUTILS_EXPORT Node curl(const Node &pSrc) const
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)
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
virtual NEKMESHUTILS_EXPORT void ProcessVertices()
Extract element vertices.
virtual NEKMESHUTILS_EXPORT void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
boost::shared_ptr< Element > ElementSharedPtr
std::pair< ModuleType, std::string > ModuleKey
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
virtual NEKMESHUTILS_EXPORT void ProcessComposites()
Generate composites.
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.