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 <<
"InputCCM: 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, triface2, triface3;
676 int quadface0, quadface1, quadface2;
680 triface0 = triface1 = triface2 = triface3 = -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)
694 else if (triface2 == -1)
699 else if (triface3 == -1)
705 if (FaceNodes[ElementFaces[i]].size() == 4)
711 else if (quadface1 == -1)
715 else if (quadface2 == -1)
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");
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");
743 int indx0, indx1, indx2, indx3, indx4;
745 indx0 = indx1 = indx2 = indx3 = indx4 = -1;
750 for (i = 0; i < 4; ++i)
752 for (j = 0; j < 3; ++j)
754 if (FaceNodes[ElementFaces[triface0]][j] ==
755 FaceNodes[ElementFaces[quadface0]][i])
765 indx2 = FaceNodes[ElementFaces[quadface0]][i];
767 else if (indx3 == -1)
769 indx3 = FaceNodes[ElementFaces[quadface0]][i];
775 "More than two vertices do not match triangular face");
782 indx0 = FaceNodes[ElementFaces[quadface0]][i];
786 indx1 = FaceNodes[ElementFaces[quadface0]][i];
792 for (
int i = 0; i < 3; ++i)
794 if ((FaceNodes[ElementFaces[triface0]][i] != indx0) &&
795 (FaceNodes[ElementFaces[triface0]][i] != indx1) &&
796 (FaceNodes[ElementFaces[triface0]][i] != indx2))
798 indx4 = FaceNodes[ElementFaces[triface0]][i];
804 Node a = *(Vnodes[indx1]) - *(Vnodes[indx0]);
806 Node b = *(Vnodes[indx4]) - *(Vnodes[indx0]);
808 Node c = *(Vnodes[indx2]) - *(Vnodes[indx0]);
812 if (acurlb_dotc < 0.0)
814 returnval[0] = indx0;
815 returnval[1] = indx1;
816 returnval[4] = indx4;
820 returnval[0] = indx1;
821 returnval[1] = indx0;
822 returnval[4] = indx4;
830 for (
int i = 0; i < 4; ++i)
832 if ((FaceNodes[ElementFaces[quadface1]][i] == returnval[1]) ||
833 (FaceNodes[ElementFaces[quadface1]][i] == indx2))
841 returnval[2] = indx2;
842 returnval[3] = indx3;
847 for (
int i = 0; i < 4; ++i)
849 if ((FaceNodes[ElementFaces[quadface2]][i] == returnval[1]) ||
850 (FaceNodes[ElementFaces[quadface2]][i] == indx2))
859 returnval[2] = indx3;
860 returnval[3] = indx2;
864 returnval[2] = indx2;
865 returnval[3] = indx3;
872 for (
int i = 0; i < 3; ++i)
874 if ((FaceNodes[ElementFaces[triface1]][i] != indx2) &&
875 (FaceNodes[ElementFaces[triface1]][i] != indx3) &&
876 (FaceNodes[ElementFaces[triface1]][i] != indx3))
878 returnval[5] = FaceNodes[ElementFaces[triface1]][i];
886 ASSERTL0(
false,
"SortFaceNodes not set up for this number of faces");
900 string fname =
m_config[
"infile"].as<
string>();
901 m_ccmErr = CCMIOOpenFile(NULL, fname.c_str(), kCCMIORead, &root);
903 CCMIOSize_t i = CCMIOSIZEC(0);
904 CCMIOID state, problem;
910 CCMIOGetState(&m_ccmErr, root, kDefaultState, &problem, &state);
911 if (m_ccmErr != kCCMIONoErr)
913 cout <<
"No state named '" << kDefaultState <<
"'" << endl;
920 CCMIONextEntity(&m_ccmErr, state, kCCMIOProcessor, &i, &
m_ccmProcessor);
925 CCMIOID mapID, vertices;
926 CCMIOSize_t nVertices, size, dims = CCMIOSIZEC(1);
930 CCMIOEntitySize(&
m_ccmErr, vertices, &nVertices, NULL);
938 int mapData[nVertices.getValue()];
939 float verts[3 * nVertices.getValue()];
940 for (
int k = 0; k < nVertices; ++k)
942 verts[3 * k] = verts[3 * k + 1] = verts[3 * k + 2] = 0.0;
952 CCMIOINDEXC(0 + nVertices));
954 &
m_ccmErr, mapID, mapData, CCMIOINDEXC(0), CCMIOINDEXC(0 + nVertices));
956 for (
int i = 0; i < nVertices; ++i)
958 Nodes.push_back(boost::shared_ptr<Node>(
959 new Node(i, verts[3 * i], verts[3 * i + 1], verts[3 * i + 2])));
968 CCMIOSize_t nFaces, size;
969 vector<int> faces, faceCells, mapData;
973 CCMIOEntitySize(&
m_ccmErr,
id, &nFaces, NULL);
975 int nf = TOINT64(nFaces);
977 faceCells.resize(2 * nf);
985 CCMIOINDEXC(kCCMIOStart),
986 CCMIOINDEXC(kCCMIOEnd));
987 faces.resize(TOINT64(size));
994 CCMIOINDEXC(kCCMIOStart),
995 CCMIOINDEXC(kCCMIOEnd));
1000 CCMIOINDEXC(kCCMIOStart),
1001 CCMIOINDEXC(kCCMIOEnd));
1005 CCMIOINDEXC(kCCMIOStart),
1006 CCMIOINDEXC(kCCMIOEnd));
1010 for (
int i = 0; i < nf; ++i)
1014 if (cnt < faces.size())
1016 int nv = faces[cnt];
1018 "Can only handle meshes with "
1019 "up to four nodes per face");
1021 for (j = 0; j < nv; ++j)
1023 if (cnt + 1 + j < faces.size())
1025 Fnodes.push_back(faces[cnt + 1 + j] - 1);
1030 FacesNodes[mapData[i] - 1] = Fnodes;
1035 for (
int i = 0; i < faceCells.size(); ++i)
1037 nelmt = max(nelmt, faceCells[i]);
1041 for (
int i = 0; i < nf; ++i)
1044 if (faceCells[2 * i])
1046 ElementFaces[faceCells[2 * i] - 1].push_back(mapData[i] - 1);
1050 if (faceCells[2 * i + 1])
1052 ElementFaces[faceCells[2 * i + 1] - 1].push_back(mapData[i] - 1);
1058 map<
int, vector<int> > &FacesNodes,
1059 Array<
OneD, vector<int> > &ElementFaces,
1060 vector<string> &Facelabels)
1063 CCMIOSize_t index = CCMIOSIZEC(0);
1065 CCMIOSize_t nFaces, size;
1066 vector<int> faces, faceCells, mapData;
1067 vector<string> facelabel;
1069 while (CCMIONextEntity(
1075 CCMIOEntitySize(&
m_ccmErr,
id, &nFaces, NULL);
1076 int nf = TOINT64(nFaces);
1078 faceCells.resize(nf);
1081 kCCMIOBoundaryFaces,
1085 CCMIOINDEXC(kCCMIOStart),
1086 CCMIOINDEXC(kCCMIOEnd));
1088 faces.resize(TOINT64(size));
1091 kCCMIOBoundaryFaces,
1095 CCMIOINDEXC(kCCMIOStart),
1096 CCMIOINDEXC(kCCMIOEnd));
1099 kCCMIOBoundaryFaces,
1101 CCMIOINDEXC(kCCMIOStart),
1102 CCMIOINDEXC(kCCMIOEnd));
1106 CCMIOINDEXC(kCCMIOStart),
1107 CCMIOINDEXC(kCCMIOEnd));
1109 CCMIOGetEntityIndex(&
m_ccmErr,
id, &boundaryVal);
1114 if (CCMIOReadOptstr(NULL,
id,
"Label", &size, NULL) == kCCMIONoErr)
1116 name =
new char[size + 1];
1117 CCMIOReadOptstr(NULL,
id,
"Label", NULL, name);
1118 Facelabels.push_back(
string(name));
1122 Facelabels.push_back(
"Not known");
1127 for (
int i = 0; i < nf; ++i)
1131 if (cnt < faces.size())
1133 int nv = faces[cnt];
1135 "Can only handle meshes with "
1136 "up to four nodes per face");
1138 for (j = 0; j < nv; ++j)
1140 if (cnt + 1 + j < faces.size())
1142 Fnodes.push_back(faces[cnt + 1 + j] - 1);
1147 FacesNodes[mapData[i] - 1] = Fnodes;
1150 vector<int> BndFaces;
1151 for (
int i = 0; i < nf; ++i)
1155 ElementFaces[faceCells[i] - 1].push_back(mapData[i] - 1);
1157 BndFaces.push_back(mapData[i] - 1);
1159 BndElementFaces.push_back(BndFaces);
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.
Represents a command-line configuration option.
NEKMESHUTILS_EXPORT Node curl(const Node &pSrc) const
std::map< std::string, ConfigOption > m_config
List of configuration values.
string msg
print "Adding",units.name, hash(units), units.description(), print "(was",id(_u),"now",id(units),")" Ensure referenced units exist
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.
static char const kDefaultState[]
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.