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;
132 nnodes = nfaces = nelements = 0;
139 boost::to_upper(line);
142 if (sscanf(line.c_str(),
"%lf", &value) == 1)
148 if ((line.find(
"NODES") != string::npos) &&
149 (line.find(
"TOTALNUMFACENODES") == string::npos))
155 start = tag.find(
"NODES=");
156 end = tag.find_first_of(
',', start);
157 nnodes = atoi(tag.substr(start + 6, end).c_str());
160 if ((line.find(
"FACES") != string::npos) &&
161 (line.find(
"NUMCONNECTEDBOUNDARYFACES") == string::npos))
167 start = tag.find(
"FACES=");
168 end = tag.find_first_of(
',', start);
169 nfaces = atoi(tag.substr(start + 6, end).c_str());
172 if (line.find(
"ELEMENTS") != string::npos)
178 start = tag.find(
"ELEMENTS=");
179 end = tag.find_first_of(
',', start);
180 nelements = atoi(tag.substr(start + 9, end).c_str());
183 if (line.find(
"ZONETYPE") != string::npos)
188 if ((line.find(
"FEPOLYGON") == string::npos) &&
189 (line.find(
"FEPOLYHEDRON") == string::npos))
192 "Routine only set up for FEPolygon or FEPolyhedron");
201 cout <<
"Setting up zone " << zcnt++;
203 vector<NekDouble> x, y, z;
206 for (i = 0; i < nnodes; ++i)
212 for (i = 0; i < nnodes; ++i)
218 for (i = 0; i < nnodes; ++i)
224 std::vector<NodeSharedPtr> Nodes;
225 for (i = 0; i < nnodes; ++i)
227 Nodes.push_back(boost::shared_ptr<Node>(
new Node(i, x[i], y[i], z[i])));
232 if (line.find(
"node count per face") == string::npos)
234 if (line.find(
"face nodes") == string::npos)
243 vector<int> Nodes_per_face;
244 if (line.find(
"node count per face") != string::npos)
247 for (i = 0; i < nfaces; ++i)
251 "Can only handle meshes with "
252 "up to four nodes per face");
253 Nodes_per_face.push_back(nodes);
260 if (line.find(
"face nodes") == string::npos)
267 vector<vector<int> > FaceNodes;
269 if (line.find(
"face nodes") != string::npos)
272 for (i = 0; i < nfaces; ++i)
276 int nodes = (Nodes_per_face.size()) ? Nodes_per_face[i] : 2;
280 for (
int j = 0; j < nodes; ++j)
285 Fnodes.push_back(nodeID - 1);
288 FaceNodes.push_back(Fnodes);
293 ASSERTL0(
false,
"Failed to find face node section");
301 if (line.find(
"left elements") == string::npos)
306 if (line.find(
"left elements") != string::npos)
310 for (i = 0; i < nfaces; ++i)
316 ElementFaces[elmtID - 1].push_back(i);
322 ASSERTL0(
false,
"Left element not found");
327 if (line.find(
"right elements") == string::npos)
332 if (line.find(
"right elements") != string::npos)
337 for (i = 0; i < nfaces; ++i)
343 ElementFaces[elmtID - 1].push_back(i);
352 ASSERTL0(
false,
"Left element not found");
355 if (Nodes_per_face.size())
357 cout <<
" (3D) " << endl;
366 for (i = 0; i < nelements; ++i)
368 if (ElementFaces[i].size() > 4)
371 Nodes, i, ElementFaces[i], FaceNodes, nComposite,
true);
378 for (i = 0; i < nelements; ++i)
380 if (ElementFaces[i].size() == 4)
383 Nodes, i, ElementFaces[i], FaceNodes, nComposite,
true);
392 cout <<
" (2D)" << endl;
396 for (i = 0; i < Nodes.size(); ++i)
400 if (it ==
m_mesh->m_vertexSet.end())
402 ASSERTL0(
false,
"Failed to find face vertex in 3D list");
410 for (i = 0; i < nelements; ++i)
412 GenElement2D(Nodes, i, ElementFaces[i], FaceNodes, nComposite);
420 map<int, int> &facelist,
421 vector<vector<int> > &FacesToPrisms,
422 vector<vector<int> > &PrismsToFaces,
423 vector<bool> &PrismDone);
427 vector<vector<int> > &FaceNodes)
431 int face1_map[3] = {0, 1, 4};
432 int face3_map[3] = {3, 2, 5};
434 map<int, bool> FacesRenumbered;
437 vector<vector<int> > FaceToPrisms(FaceNodes.size());
438 vector<vector<int> > PrismToFaces(ElementFaces.num_elements());
439 map<int, int> Prisms;
444 for (i = 0; i < ElementFaces.num_elements(); ++i)
447 if (ElementFaces[i].size() == 5)
449 vector<int> LocTriFaces;
451 for (j = 0; j < ElementFaces[i].size(); ++j)
453 if (FaceNodes[ElementFaces[i][j]].size() == 3)
455 LocTriFaces.push_back(j);
459 if (LocTriFaces.size() == 2)
463 PrismToFaces[i].push_back(ElementFaces[i][LocTriFaces[0]]);
464 PrismToFaces[i].push_back(ElementFaces[i][LocTriFaces[1]]);
466 FaceToPrisms[ElementFaces[i][LocTriFaces[0]]].push_back(i);
467 FaceToPrisms[ElementFaces[i][LocTriFaces[1]]].push_back(i);
472 vector<bool> FacesDone(FaceNodes.size(),
false);
473 vector<bool> PrismDone(ElementFaces.num_elements(),
false);
478 for (PrismIt = Prisms.begin(); PrismIt != Prisms.end(); ++PrismIt)
480 int elmtid = PrismIt->first;
481 map<int, int> facelist;
484 if (PrismDone[elmtid])
492 elmtid, facelist, FaceToPrisms, PrismToFaces, PrismDone);
495 for (faceIt = facelist.begin(); faceIt != facelist.end(); faceIt++)
497 int faceid = faceIt->second;
499 for (i = 0; i < FaceToPrisms[faceid].size(); ++i)
501 int prismid = FaceToPrisms[faceid][i];
503 if ((FacesDone[PrismToFaces[prismid][0]] ==
true) &&
504 (FacesDone[PrismToFaces[prismid][1]] ==
true))
512 if ((FacesDone[PrismToFaces[prismid][0]] ==
false) &&
513 (FacesDone[PrismToFaces[prismid][1]] ==
false))
517 for (i = 0; i < 3; ++i)
519 if (NodeReordering[Nodes[face1_map[i]]] == -1)
521 NodeReordering[Nodes[face1_map[i]]] = nodeid++;
525 for (i = 0; i < 3; ++i)
527 if (NodeReordering[Nodes[face3_map[i]]] == -1)
529 NodeReordering[Nodes[face3_map[i]]] = nodeid++;
533 else if ((FacesDone[PrismToFaces[prismid][0]] ==
false) &&
534 (FacesDone[PrismToFaces[prismid][1]] ==
true))
537 int max_id1, max_id2;
539 max_id1 = (NodeReordering[Nodes[face3_map[0]]] <
540 NodeReordering[Nodes[face3_map[1]]])
543 max_id2 = (NodeReordering[Nodes[face3_map[max_id1]]] <
544 NodeReordering[Nodes[face3_map[2]]])
549 int id0 = (max_id1 == 1) ? 0 : 1;
551 if (NodeReordering[Nodes[face1_map[id0]]] == -1)
553 NodeReordering[Nodes[face1_map[id0]]] = nodeid++;
556 if (NodeReordering[Nodes[face1_map[max_id1]]] == -1)
558 NodeReordering[Nodes[face1_map[max_id1]]] =
562 if (NodeReordering[Nodes[face1_map[max_id2]]] == -1)
564 NodeReordering[Nodes[face1_map[max_id2]]] =
568 else if ((FacesDone[PrismToFaces[prismid][0]] ==
true) &&
569 (FacesDone[PrismToFaces[prismid][1]] ==
false))
572 int max_id1, max_id2;
574 max_id1 = (NodeReordering[Nodes[face1_map[0]]] <
575 NodeReordering[Nodes[face1_map[1]]])
578 max_id2 = (NodeReordering[Nodes[face1_map[max_id1]]] <
579 NodeReordering[Nodes[face1_map[2]]])
584 int id0 = (max_id1 == 1) ? 0 : 1;
586 if (NodeReordering[Nodes[face3_map[id0]]] == -1)
588 NodeReordering[Nodes[face3_map[id0]]] = nodeid++;
591 if (NodeReordering[Nodes[face3_map[max_id1]]] == -1)
593 NodeReordering[Nodes[face3_map[max_id1]]] =
597 if (NodeReordering[Nodes[face3_map[max_id2]]] == -1)
599 NodeReordering[Nodes[face3_map[max_id2]]] =
609 for (i = 0; i < NodeReordering.num_elements(); ++i)
611 if (NodeReordering[i] == -1)
613 NodeReordering[i] = nodeid++;
617 ASSERTL1(nodeid == NodeReordering.num_elements(),
618 "Have not renumbered all nodes");
621 for (i = 0; i < FaceNodes.size(); ++i)
623 for (j = 0; j < FaceNodes[i].size(); ++j)
625 FaceNodes[i][j] = NodeReordering[FaceNodes[i][j]];
629 vector<NodeSharedPtr> save(Vnodes);
630 for (i = 0; i < Vnodes.size(); ++i)
632 Vnodes[NodeReordering[i]] = save[i];
633 Vnodes[NodeReordering[i]]->SetID(NodeReordering[i]);
638 map<int, int> &facelist,
639 vector<vector<int> > &FaceToPrisms,
640 vector<vector<int> > &PrismToFaces,
641 vector<bool> &PrismDone)
643 if (PrismDone[prismid] ==
false)
645 PrismDone[prismid] =
true;
648 int face = PrismToFaces[prismid][0];
649 facelist[face] = face;
650 for (
int i = 0; i < FaceToPrisms[face].size(); ++i)
660 face = PrismToFaces[prismid][1];
661 facelist[face] = face;
662 for (
int i = 0; i < FaceToPrisms[face].size(); ++i)
675 vector<int> &ElementFaces,
676 vector<vector<int> > &FaceNodes,
682 if (ElementFaces.size() == 3)
686 else if (ElementFaces.size() == 4)
692 ASSERTL0(
false,
"Not set up for elements which are not Tets or Prism");
697 tags.push_back(nComposite);
700 vector<NodeSharedPtr> nodeList;
702 for (
int j = 0; j < Nodes.num_elements(); ++j)
704 nodeList.push_back(VertNodes[Nodes[j]]);
712 m_mesh->m_element[E->GetDim()].push_back(E);
717 vector<int> &ElementFaces,
718 vector<vector<int> > &FaceNodes,
725 int nnodes = Nodes.num_elements();
726 map<LibUtilities::ShapeType, int> domainComposite;
736 else if (nnodes == 5)
740 else if (nnodes == 6)
747 ASSERTL0(
false,
"Not set up for elements which are not Tets or Prism");
752 tags.push_back(nComposite);
755 vector<NodeSharedPtr> nodeList;
756 for (
int j = 0; j < Nodes.num_elements(); ++j)
758 nodeList.push_back(VertNodes[Nodes[j]]);
764 ElmtConfig conf(elType, 1,
true,
true, DoOrient);
768 m_mesh->m_element[E->GetDim()].push_back(E);
772 cout <<
"Warning: Pyramid detected " << endl;
777 vector<int> &ElementFaces,
778 vector<vector<int> > &FaceNodes)
783 if (ElementFaces.size() == 3)
787 returnval[0] = FaceNodes[ElementFaces[0]][0];
788 returnval[1] = FaceNodes[ElementFaces[0]][1];
791 for (i = 0; i < 2; ++i)
793 if ((FaceNodes[ElementFaces[1]][i] != returnval[0]) &&
794 (FaceNodes[ElementFaces[1]][i] != returnval[1]))
796 returnval[2] = FaceNodes[ElementFaces[1]][i];
801 else if (ElementFaces.size() == 4)
805 int indx0 = FaceNodes[ElementFaces[0]][0];
806 int indx1 = FaceNodes[ElementFaces[0]][1];
811 for (j = 1; j < 4; ++j)
813 for (i = 0; i < 2; ++i)
815 if ((FaceNodes[ElementFaces[j]][i] != indx0) &&
816 (FaceNodes[ElementFaces[j]][i] != indx1))
820 indx2 = FaceNodes[ElementFaces[j]][i];
822 else if (indx2 != -1)
824 if (FaceNodes[ElementFaces[j]][i] != indx2)
826 indx3 = FaceNodes[ElementFaces[j]][i];
833 ASSERTL1((indx2 != -1) && (indx3 != -1),
834 "Failed to find vertex 3 or 4");
837 Node a = *(Vnodes[indx1]) - *(Vnodes[indx0]);
839 Node b = *(Vnodes[indx2]) - *(Vnodes[indx0]);
843 Node c = *(Vnodes[indx1]) - *(Vnodes[indx2]);
845 Node d = *(Vnodes[indx3]) - *(Vnodes[indx2]);
849 if (acurlb_dot_acurld > 0.0)
851 returnval[0] = indx0;
852 returnval[1] = indx1;
853 returnval[2] = indx2;
854 returnval[3] = indx3;
858 returnval[0] = indx0;
859 returnval[1] = indx1;
860 returnval[2] = indx3;
861 returnval[3] = indx2;
869 vector<int> &ElementFaces,
870 vector<vector<int> > &FaceNodes)
876 if (ElementFaces.size() == 4)
878 ASSERTL1(FaceNodes[ElementFaces[0]].size() == 3,
879 "Face is not triangular");
883 int indx0 = FaceNodes[ElementFaces[0]][0];
884 int indx1 = FaceNodes[ElementFaces[0]][1];
885 int indx2 = FaceNodes[ElementFaces[0]][2];
889 Node a = *(Vnodes[indx1]) - *(Vnodes[indx0]);
891 Node b = *(Vnodes[indx2]) - *(Vnodes[indx0]);
894 ASSERTL1(FaceNodes[ElementFaces[1]].size() == 3,
895 "Face is not triangular");
896 for (i = 0; i < 3; ++i)
899 if ((FaceNodes[ElementFaces[1]][i] != indx0) &&
900 (FaceNodes[ElementFaces[1]][i] != indx1) &&
901 (FaceNodes[ElementFaces[1]][i] != indx2))
903 indx3 = FaceNodes[ElementFaces[1]][i];
909 Node c = *(Vnodes[indx3]) - *(Vnodes[indx0]);
913 if (acurlb_dotc < 0.0)
915 returnval[0] = indx0;
916 returnval[1] = indx1;
917 returnval[2] = indx2;
918 returnval[3] = indx3;
922 returnval[0] = indx1;
923 returnval[1] = indx0;
924 returnval[2] = indx2;
925 returnval[3] = indx3;
928 else if (ElementFaces.size() == 5)
930 int triface0, triface1;
931 int quadface0, quadface1, quadface2;
935 triface0 = triface1 = -1;
936 quadface0 = quadface1 = quadface2 = -1;
937 for (i = 0; i < 5; ++i)
939 if (FaceNodes[ElementFaces[i]].size() == 3)
945 else if (triface1 == -1)
955 if (FaceNodes[ElementFaces[i]].size() == 4)
961 else if (quadface1 == -1)
965 else if (quadface2 == -1)
982 int indx0, indx1, indx2, indx3, indx4;
984 indx0 = indx1 = indx2 = indx3 = indx4 = -1;
989 for (i = 0; i < 4; ++i)
991 for (j = 0; j < 3; ++j)
993 if (FaceNodes[ElementFaces[triface0]][j] ==
994 FaceNodes[ElementFaces[quadface0]][i])
1004 indx2 = FaceNodes[ElementFaces[quadface0]][i];
1006 else if (indx3 == -1)
1008 indx3 = FaceNodes[ElementFaces[quadface0]][i];
1014 "More than two vertices do not match triangular face");
1021 indx0 = FaceNodes[ElementFaces[quadface0]][i];
1025 indx1 = FaceNodes[ElementFaces[quadface0]][i];
1031 for (
int i = 0; i < 3; ++i)
1033 if ((FaceNodes[ElementFaces[triface0]][i] != indx0) &&
1034 (FaceNodes[ElementFaces[triface0]][i] != indx1) &&
1035 (FaceNodes[ElementFaces[triface0]][i] != indx2))
1037 indx4 = FaceNodes[ElementFaces[triface0]][i];
1043 Node a = *(Vnodes[indx1]) - *(Vnodes[indx0]);
1045 Node b = *(Vnodes[indx4]) - *(Vnodes[indx0]);
1047 Node c = *(Vnodes[indx2]) - *(Vnodes[indx0]);
1051 if (acurlb_dotc < 0.0)
1053 returnval[0] = indx0;
1054 returnval[1] = indx1;
1055 returnval[4] = indx4;
1059 returnval[0] = indx1;
1060 returnval[1] = indx0;
1061 returnval[4] = indx4;
1069 for (
int i = 0; i < 4; ++i)
1071 if ((FaceNodes[ElementFaces[quadface1]][i] == returnval[1]) ||
1072 (FaceNodes[ElementFaces[quadface1]][i] == indx2))
1080 returnval[2] = indx2;
1081 returnval[3] = indx3;
1086 for (
int i = 0; i < 4; ++i)
1088 if ((FaceNodes[ElementFaces[quadface2]][i] == returnval[1]) ||
1089 (FaceNodes[ElementFaces[quadface2]][i] == indx2))
1098 returnval[2] = indx3;
1099 returnval[3] = indx2;
1103 returnval[2] = indx2;
1104 returnval[3] = indx3;
1108 if (isPrism ==
true)
1111 for (
int i = 0; i < 3; ++i)
1113 if ((FaceNodes[ElementFaces[triface1]][i] != indx2) &&
1114 (FaceNodes[ElementFaces[triface1]][i] != indx3) &&
1115 (FaceNodes[ElementFaces[triface1]][i] != indx3))
1117 returnval[5] = FaceNodes[ElementFaces[triface1]][i];
1125 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.
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.
MeshSharedPtr m_mesh
Mesh object.
ElementFactory & GetElementFactory()
Represents a point in the domain.
virtual void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
virtual void ProcessVertices()
Extract element vertices.
virtual void ProcessElements()
Generate element IDs.
NEKMESHUTILS_EXPORT Node curl(const Node &pSrc) const
virtual void ProcessComposites()
Generate composites.
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 void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
boost::shared_ptr< Element > ElementSharedPtr
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.