36 #include <boost/algorithm/string.hpp>
54 "Reads mesh from Star CCM (.ccm).");
61 "Just write out tags from star file for each surface/composite");
84 cout <<
"InputStarTec: Start reading file..." << endl;
108 std::vector<NodeSharedPtr> Nodes;
112 map<int, vector<int> > FaceNodes;
119 vector<vector<int> > BndElementFaces;
120 vector<string> Facelabels;
123 if (
m_config[
"writelabelsonly"].beenSet)
127 for (i = 0; i < BndElementFaces.size(); ++i)
129 cout <<
" 2D Zone (composite = " << nComposite
130 <<
", label = " << Facelabels[i] <<
")" << endl;
144 int nelements = ElementFaces.num_elements();
145 cout <<
" Generating 3D Zones: ";
147 for (i = 0; i < nelements; ++i)
150 if (ElementFaces[i].size() > 4)
153 Nodes, i, ElementFaces[i], FaceNodes, nComposite,
true);
157 cout << cnt <<
" Prisms,";
163 for (i = 0; i < nelements; ++i)
165 if (ElementFaces[i].size() == 4)
168 Nodes, i, ElementFaces[i], FaceNodes, nComposite,
true);
172 cout << cnt <<
" Tets" << endl;
178 for (i = 0; i < BndElementFaces.size(); ++i)
180 cout <<
" Generating 2D Zone (composite = " << nComposite
181 <<
", label = " << Facelabels[i] <<
")" << endl;
183 for (
int j = 0; j < BndElementFaces[i].size(); ++j)
186 if (FaceNodes.count(BndElementFaces[i][j]))
189 Nodes, j, FaceNodes[BndElementFaces[i][j]], nComposite);
193 string msg =
"Failed to find FaceNodes for Face ";
194 msg += boost::lexical_cast<
string>(BndElementFaces[i][j]);
199 m_mesh->m_faceLabels[nComposite] = Facelabels[i];
205 map<int, int> &facelist,
206 vector<vector<int> > &FacesToPrisms,
207 vector<vector<int> > &PrismsToFaces,
208 vector<bool> &PrismDone);
212 map<
int, vector<int> > &FaceNodes)
216 int face1_map[3] = {0, 1, 4};
217 int face3_map[3] = {3, 2, 5};
219 map<int, bool> FacesRenumbered;
222 vector<vector<int> > FaceToPrisms(FaceNodes.size());
223 vector<vector<int> > PrismToFaces(ElementFaces.num_elements());
224 map<int, int> Prisms;
229 for (i = 0; i < ElementFaces.num_elements(); ++i)
232 if (ElementFaces[i].size() == 5)
234 vector<int> LocTriFaces;
236 for (j = 0; j < ElementFaces[i].size(); ++j)
238 if (FaceNodes[ElementFaces[i][j]].size() == 3)
240 LocTriFaces.push_back(j);
244 if (LocTriFaces.size() == 2)
248 PrismToFaces[i].push_back(ElementFaces[i][LocTriFaces[0]]);
249 PrismToFaces[i].push_back(ElementFaces[i][LocTriFaces[1]]);
251 FaceToPrisms[ElementFaces[i][LocTriFaces[0]]].push_back(i);
252 FaceToPrisms[ElementFaces[i][LocTriFaces[1]]].push_back(i);
257 vector<bool> FacesDone(FaceNodes.size(),
false);
258 vector<bool> PrismDone(ElementFaces.num_elements(),
false);
263 for (PrismIt = Prisms.begin(); PrismIt != Prisms.end(); ++PrismIt)
265 int elmtid = PrismIt->first;
266 map<int, int> facelist;
269 if (PrismDone[elmtid])
277 elmtid, facelist, FaceToPrisms, PrismToFaces, PrismDone);
280 for (faceIt = facelist.begin(); faceIt != facelist.end(); faceIt++)
282 int faceid = faceIt->second;
284 for (i = 0; i < FaceToPrisms[faceid].size(); ++i)
286 int prismid = FaceToPrisms[faceid][i];
288 if ((FacesDone[PrismToFaces[prismid][0]] ==
true) &&
289 (FacesDone[PrismToFaces[prismid][1]] ==
true))
297 if ((FacesDone[PrismToFaces[prismid][0]] ==
false) &&
298 (FacesDone[PrismToFaces[prismid][1]] ==
false))
302 for (i = 0; i < 3; ++i)
304 if (NodeReordering[Nodes[face1_map[i]]] == -1)
306 NodeReordering[Nodes[face1_map[i]]] = nodeid++;
310 for (i = 0; i < 3; ++i)
312 if (NodeReordering[Nodes[face3_map[i]]] == -1)
314 NodeReordering[Nodes[face3_map[i]]] = nodeid++;
318 else if ((FacesDone[PrismToFaces[prismid][0]] ==
false) &&
319 (FacesDone[PrismToFaces[prismid][1]] ==
true))
322 int max_id1, max_id2;
324 max_id1 = (NodeReordering[Nodes[face3_map[0]]] <
325 NodeReordering[Nodes[face3_map[1]]])
328 max_id2 = (NodeReordering[Nodes[face3_map[max_id1]]] <
329 NodeReordering[Nodes[face3_map[2]]])
334 int id0 = (max_id1 == 1) ? 0 : 1;
336 if (NodeReordering[Nodes[face1_map[id0]]] == -1)
338 NodeReordering[Nodes[face1_map[id0]]] = nodeid++;
341 if (NodeReordering[Nodes[face1_map[max_id1]]] == -1)
343 NodeReordering[Nodes[face1_map[max_id1]]] =
347 if (NodeReordering[Nodes[face1_map[max_id2]]] == -1)
349 NodeReordering[Nodes[face1_map[max_id2]]] =
353 else if ((FacesDone[PrismToFaces[prismid][0]] ==
true) &&
354 (FacesDone[PrismToFaces[prismid][1]] ==
false))
357 int max_id1, max_id2;
359 max_id1 = (NodeReordering[Nodes[face1_map[0]]] <
360 NodeReordering[Nodes[face1_map[1]]])
363 max_id2 = (NodeReordering[Nodes[face1_map[max_id1]]] <
364 NodeReordering[Nodes[face1_map[2]]])
369 int id0 = (max_id1 == 1) ? 0 : 1;
371 if (NodeReordering[Nodes[face3_map[id0]]] == -1)
373 NodeReordering[Nodes[face3_map[id0]]] = nodeid++;
376 if (NodeReordering[Nodes[face3_map[max_id1]]] == -1)
378 NodeReordering[Nodes[face3_map[max_id1]]] =
382 if (NodeReordering[Nodes[face3_map[max_id2]]] == -1)
384 NodeReordering[Nodes[face3_map[max_id2]]] =
394 for (i = 0; i < NodeReordering.num_elements(); ++i)
396 if (NodeReordering[i] == -1)
398 NodeReordering[i] = nodeid++;
402 ASSERTL1(nodeid == NodeReordering.num_elements(),
403 "Have not renumbered all nodes");
406 for (i = 0; i < FaceNodes.size(); ++i)
408 for (j = 0; j < FaceNodes[i].size(); ++j)
410 FaceNodes[i][j] = NodeReordering[FaceNodes[i][j]];
414 vector<NodeSharedPtr> save(Vnodes);
415 for (i = 0; i < Vnodes.size(); ++i)
417 Vnodes[NodeReordering[i]] = save[i];
418 Vnodes[NodeReordering[i]]->SetID(NodeReordering[i]);
423 map<int, int> &facelist,
424 vector<vector<int> > &FaceToPrisms,
425 vector<vector<int> > &PrismToFaces,
426 vector<bool> &PrismDone)
428 if (PrismDone[prismid] ==
false)
430 PrismDone[prismid] =
true;
433 int face = PrismToFaces[prismid][0];
434 facelist[face] = face;
435 for (
int i = 0; i < FaceToPrisms[face].size(); ++i)
445 face = PrismToFaces[prismid][1];
446 facelist[face] = face;
447 for (
int i = 0; i < FaceToPrisms[face].size(); ++i)
460 vector<int> &FaceNodes,
465 if (FaceNodes.size() == 3)
469 else if (FaceNodes.size() == 4)
475 ASSERTL0(
false,
"Not set up for elements which are not Tets or Prism");
480 tags.push_back(nComposite);
483 vector<NodeSharedPtr> nodeList;
485 for (
int j = 0; j < Nodes.num_elements(); ++j)
487 nodeList.push_back(VertNodes[Nodes[j]]);
495 m_mesh->m_element[E->GetDim()].push_back(E);
500 vector<int> &ElementFaces,
501 map<
int, vector<int> > &FaceNodes,
508 int nnodes = Nodes.num_elements();
509 map<LibUtilities::ShapeType, int> domainComposite;
519 else if (nnodes == 5)
523 else if (nnodes == 6)
530 ASSERTL0(
false,
"Not set up for elements which are not Tets or Prism");
535 tags.push_back(nComposite);
538 vector<NodeSharedPtr> nodeList;
539 for (
int j = 0; j < Nodes.num_elements(); ++j)
541 nodeList.push_back(VertNodes[Nodes[j]]);
547 ElmtConfig conf(elType, 1,
true,
true, DoOrient);
551 m_mesh->m_element[E->GetDim()].push_back(E);
555 cout <<
"Warning: Pyramid detected " << endl;
560 vector<int> &FaceNodes)
564 if (FaceNodes.size() == 3)
568 returnval[0] = FaceNodes[0];
569 returnval[1] = FaceNodes[1];
570 returnval[2] = FaceNodes[2];
572 else if (FaceNodes.size() == 4)
576 int indx0 = FaceNodes[0];
577 int indx1 = FaceNodes[1];
578 int indx2 = FaceNodes[2];
579 int indx3 = FaceNodes[3];
582 Node a = *(Vnodes[indx1]) - *(Vnodes[indx0]);
584 Node b = *(Vnodes[indx2]) - *(Vnodes[indx0]);
588 Node c = *(Vnodes[indx1]) - *(Vnodes[indx2]);
590 Node d = *(Vnodes[indx3]) - *(Vnodes[indx2]);
594 if (acurlb_dot_acurld > 0.0)
596 returnval[0] = indx0;
597 returnval[1] = indx1;
598 returnval[2] = indx2;
599 returnval[3] = indx3;
603 returnval[0] = indx0;
604 returnval[1] = indx1;
605 returnval[2] = indx3;
606 returnval[3] = indx2;
614 vector<int> &ElementFaces,
615 map<
int, vector<int> > &FaceNodes)
621 if (ElementFaces.size() == 4)
623 ASSERTL1(FaceNodes[ElementFaces[0]].size() == 3,
624 "Face is not triangular");
628 int indx0 = FaceNodes[ElementFaces[0]][0];
629 int indx1 = FaceNodes[ElementFaces[0]][1];
630 int indx2 = FaceNodes[ElementFaces[0]][2];
634 Node a = *(Vnodes[indx1]) - *(Vnodes[indx0]);
636 Node b = *(Vnodes[indx2]) - *(Vnodes[indx0]);
639 ASSERTL1(FaceNodes[ElementFaces[1]].size() == 3,
640 "Face is not triangular");
641 for (i = 0; i < 3; ++i)
644 if ((FaceNodes[ElementFaces[1]][i] != indx0) &&
645 (FaceNodes[ElementFaces[1]][i] != indx1) &&
646 (FaceNodes[ElementFaces[1]][i] != indx2))
648 indx3 = FaceNodes[ElementFaces[1]][i];
654 Node c = *(Vnodes[indx3]) - *(Vnodes[indx0]);
658 if (acurlb_dotc < 0.0)
660 returnval[0] = indx0;
661 returnval[1] = indx1;
662 returnval[2] = indx2;
663 returnval[3] = indx3;
667 returnval[0] = indx1;
668 returnval[1] = indx0;
669 returnval[2] = indx2;
670 returnval[3] = indx3;
673 else if (ElementFaces.size() == 5)
675 int triface0, triface1;
676 int quadface0, quadface1, quadface2;
680 triface0 = triface1 = -1;
681 quadface0 = quadface1 = quadface2 = -1;
682 for (i = 0; i < 5; ++i)
684 if (FaceNodes[ElementFaces[i]].size() == 3)
690 else if (triface1 == -1)
700 if (FaceNodes[ElementFaces[i]].size() == 4)
706 else if (quadface1 == -1)
710 else if (quadface2 == -1)
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");
731 cout <<
"Pyramid found with vertices: " << endl;
732 for (i = 0; i < 5; ++i)
734 for (j = 0; j < FaceNodes[ElementFaces[i]].size(); ++j)
736 vertids.insert(FaceNodes[ElementFaces[i]][j]);
739 for (it = vertids.begin(); it != vertids.end(); ++it)
741 cout << Vnodes[*it] << endl;
744 ASSERTL0(
false,
"Not yet set up for pyramids");
749 int indx0, indx1, indx2, indx3, indx4;
751 indx0 = indx1 = indx2 = indx3 = indx4 = -1;
756 for (i = 0; i < 4; ++i)
758 for (j = 0; j < 3; ++j)
760 if (FaceNodes[ElementFaces[triface0]][j] ==
761 FaceNodes[ElementFaces[quadface0]][i])
771 indx2 = FaceNodes[ElementFaces[quadface0]][i];
773 else if (indx3 == -1)
775 indx3 = FaceNodes[ElementFaces[quadface0]][i];
781 "More than two vertices do not match triangular face");
788 indx0 = FaceNodes[ElementFaces[quadface0]][i];
792 indx1 = FaceNodes[ElementFaces[quadface0]][i];
798 for (
int i = 0; i < 3; ++i)
800 if ((FaceNodes[ElementFaces[triface0]][i] != indx0) &&
801 (FaceNodes[ElementFaces[triface0]][i] != indx1) &&
802 (FaceNodes[ElementFaces[triface0]][i] != indx2))
804 indx4 = FaceNodes[ElementFaces[triface0]][i];
810 Node a = *(Vnodes[indx1]) - *(Vnodes[indx0]);
812 Node b = *(Vnodes[indx4]) - *(Vnodes[indx0]);
814 Node c = *(Vnodes[indx2]) - *(Vnodes[indx0]);
818 if (acurlb_dotc < 0.0)
820 returnval[0] = indx0;
821 returnval[1] = indx1;
822 returnval[4] = indx4;
826 returnval[0] = indx1;
827 returnval[1] = indx0;
828 returnval[4] = indx4;
836 for (
int i = 0; i < 4; ++i)
838 if ((FaceNodes[ElementFaces[quadface1]][i] == returnval[1]) ||
839 (FaceNodes[ElementFaces[quadface1]][i] == indx2))
847 returnval[2] = indx2;
848 returnval[3] = indx3;
853 for (
int i = 0; i < 4; ++i)
855 if ((FaceNodes[ElementFaces[quadface2]][i] == returnval[1]) ||
856 (FaceNodes[ElementFaces[quadface2]][i] == indx2))
865 returnval[2] = indx3;
866 returnval[3] = indx2;
870 returnval[2] = indx2;
871 returnval[3] = indx3;
878 for (
int i = 0; i < 3; ++i)
880 if ((FaceNodes[ElementFaces[triface1]][i] != indx2) &&
881 (FaceNodes[ElementFaces[triface1]][i] != indx3) &&
882 (FaceNodes[ElementFaces[triface1]][i] != indx3))
884 returnval[5] = FaceNodes[ElementFaces[triface1]][i];
892 ASSERTL0(
false,
"SortFaceNodes not set up for this number of faces");
906 string fname =
m_config[
"infile"].as<
string>();
907 m_ccmErr = CCMIOOpenFile(NULL, fname.c_str(), kCCMIORead, &root);
909 CCMIOSize_t i = CCMIOSIZEC(0);
910 CCMIOID state, problem;
916 CCMIOGetState(&m_ccmErr, root, kDefaultState, &problem, &state);
917 if (m_ccmErr != kCCMIONoErr)
919 cout <<
"No state named '" << kDefaultState <<
"'" << endl;
926 CCMIONextEntity(&m_ccmErr, state, kCCMIOProcessor, &i, &
m_ccmProcessor);
931 CCMIOID mapID, vertices;
932 CCMIOSize_t nVertices, size, dims = CCMIOSIZEC(1);
936 CCMIOEntitySize(&
m_ccmErr, vertices, &nVertices, NULL);
944 int mapData[nVertices.getValue()];
945 float verts[3 * nVertices.getValue()];
946 for (
int k = 0; k < nVertices; ++k)
948 verts[3 * k] = verts[3 * k + 1] = verts[3 * k + 2] = 0.0;
958 CCMIOINDEXC(0 + nVertices));
960 &
m_ccmErr, mapID, mapData, CCMIOINDEXC(0), CCMIOINDEXC(0 + nVertices));
962 for (
int i = 0; i < nVertices; ++i)
964 Nodes.push_back(boost::shared_ptr<Node>(
965 new Node(i, verts[3 * i], verts[3 * i + 1], verts[3 * i + 2])));
974 CCMIOSize_t nFaces, size;
975 vector<int> faces, faceCells, mapData;
979 CCMIOEntitySize(&
m_ccmErr,
id, &nFaces, NULL);
981 int nf = TOINT64(nFaces);
983 faceCells.resize(2 * nf);
991 CCMIOINDEXC(kCCMIOStart),
992 CCMIOINDEXC(kCCMIOEnd));
993 faces.resize(TOINT64(size));
1000 CCMIOINDEXC(kCCMIOStart),
1001 CCMIOINDEXC(kCCMIOEnd));
1004 kCCMIOInternalFaces,
1006 CCMIOINDEXC(kCCMIOStart),
1007 CCMIOINDEXC(kCCMIOEnd));
1011 CCMIOINDEXC(kCCMIOStart),
1012 CCMIOINDEXC(kCCMIOEnd));
1016 for (
int i = 0; i < nf; ++i)
1020 if (cnt < faces.size())
1022 int nv = faces[cnt];
1024 "Can only handle meshes with "
1025 "up to four nodes per face");
1027 for (j = 0; j < nv; ++j)
1029 if (cnt + 1 + j < faces.size())
1031 Fnodes.push_back(faces[cnt + 1 + j] - 1);
1036 FacesNodes[mapData[i] - 1] = Fnodes;
1041 for (
int i = 0; i < faceCells.size(); ++i)
1043 nelmt = max(nelmt, faceCells[i]);
1047 for (
int i = 0; i < nf; ++i)
1050 if (faceCells[2 * i])
1052 ElementFaces[faceCells[2 * i] - 1].push_back(mapData[i] - 1);
1056 if (faceCells[2 * i + 1])
1058 ElementFaces[faceCells[2 * i + 1] - 1].push_back(mapData[i] - 1);
1064 map<
int, vector<int> > &FacesNodes,
1065 Array<
OneD, vector<int> > &ElementFaces,
1066 vector<string> &Facelabels)
1069 CCMIOSize_t index = CCMIOSIZEC(0);
1071 CCMIOSize_t nFaces, size;
1072 vector<int> faces, faceCells, mapData;
1073 vector<string> facelabel;
1075 while (CCMIONextEntity(
1081 CCMIOEntitySize(&
m_ccmErr,
id, &nFaces, NULL);
1082 int nf = TOINT64(nFaces);
1084 faceCells.resize(nf);
1087 kCCMIOBoundaryFaces,
1091 CCMIOINDEXC(kCCMIOStart),
1092 CCMIOINDEXC(kCCMIOEnd));
1094 faces.resize(TOINT64(size));
1097 kCCMIOBoundaryFaces,
1101 CCMIOINDEXC(kCCMIOStart),
1102 CCMIOINDEXC(kCCMIOEnd));
1105 kCCMIOBoundaryFaces,
1107 CCMIOINDEXC(kCCMIOStart),
1108 CCMIOINDEXC(kCCMIOEnd));
1112 CCMIOINDEXC(kCCMIOStart),
1113 CCMIOINDEXC(kCCMIOEnd));
1115 CCMIOGetEntityIndex(&
m_ccmErr,
id, &boundaryVal);
1120 if (CCMIOReadOptstr(NULL,
id,
"Label", &size, NULL) == kCCMIONoErr)
1122 name =
new char[size + 1];
1123 CCMIOReadOptstr(NULL,
id,
"Label", NULL, name);
1124 Facelabels.push_back(
string(name));
1128 Facelabels.push_back(
"Not known");
1133 for (
int i = 0; i < nf; ++i)
1137 if (cnt < faces.size())
1139 int nv = faces[cnt];
1141 "Can only handle meshes with "
1142 "up to four nodes per face");
1144 for (j = 0; j < nv; ++j)
1146 if (cnt + 1 + j < faces.size())
1148 Fnodes.push_back(faces[cnt + 1 + j] - 1);
1153 FacesNodes[mapData[i] - 1] = Fnodes;
1156 vector<int> BndFaces;
1157 for (
int i = 0; i < nf; ++i)
1161 ElementFaces[faceCells[i] - 1].push_back(mapData[i] - 1);
1163 BndFaces.push_back(mapData[i] - 1);
1165 BndElementFaces.push_back(BndFaces);
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.
map< string, ConfigOption > m_config
List of configuration values.
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.
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)
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
virtual void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
static char const kDefaultState[]
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.