70 static label MSHLINE = 1;
72 static label MSHTRI = 2;
73 static label MSHQUAD = 3;
74 static label MSHTET = 4;
77 static label MSHHEX = 5;
78 static label MSHPRISM = 6;
79 static label MSHPYR = 7;
95 while (
line.size() < 4 ||
line.substr(0, 4) !=
"$End");
120 if (!meshPointMap.found(meshF[0]))
122 Warning<<
"Not using gmsh face " << meshF
123 <<
" since zero vertex is not on boundary of polyMesh" <<
endl;
136 const face&
f = pp[facei];
143 if (meshF.found(
f[fp]))
149 if (nMatched == meshF.size())
175 if (meshF.found(
f[fp]))
181 if (nMatched == meshF.size())
200 for (
const face&
f : faces)
205 if (((
points[
f[0]] - cc) & areaNorm) < 0)
226 const auto zoneFnd = physToZone.cfind(regPhys);
231 zoneCells[zoneFnd()].append(celli);
236 const label zonei = zoneCells.size();
237 zoneCells.setSize(zonei+1);
240 Info<<
"Mapping region " << regPhys <<
" to Foam cellZone "
242 physToZone.insert(regPhys, zonei);
244 zoneToPhys[zonei] = regPhys;
245 zoneCells[zonei].
append(celli);
251 scalar readMeshFormat(
IFstream& inFile)
253 Info<<
"Starting to read mesh format at line "
262 label asciiFlag, nBytes;
263 lineStr >>
version >> asciiFlag >> nBytes;
265 Info<<
"Read format version " <<
version <<
" ascii " << asciiFlag <<
endl;
270 <<
"Can only read ascii msh files."
278 if (tag !=
"$EndMeshFormat")
281 <<
"Did not find $ENDNOD tag on line "
302 Info<<
"Vertices to be read: " << nVerts <<
endl;
306 for (label pointi = 0; pointi < nVerts; pointi++)
309 scalar xVal, yVal, zVal;
315 lineStr >> mshLabel >> xVal >> yVal >> zVal;
323 mshToFoam.insert(mshLabel, pointi);
326 Info<<
"Vertices read:" << mshToFoam.size() <<
endl;
332 if (tag !=
"$ENDNOD" && tag !=
"$EndNodes")
335 <<
"Did not find $ENDNOD tag on line "
351 label nEntities, nVerts;
352 lineStr >> nEntities >> nVerts;
354 Info<<
"Vertices to be read: " << nVerts <<
endl;
362 for (label entityi = 0; entityi < nEntities; entityi++)
364 label entityDim, entityLabel, isParametric, nNodes;
365 scalar xVal, yVal, zVal;
370 lineStr >> entityDim >> entityLabel >> isParametric >> nNodes;
374 for (label eNode = 0; eNode < nNodes; ++eNode)
378 lineStr >> nodeIDs[eNode];
382 for (label eNode = 0; eNode < nNodes; ++eNode)
386 lineStr >> xVal >> yVal >> zVal;
391 mshToFoam.insert(nodeIDs[eNode], pointi++);
396 Info<<
"Vertices read: " << mshToFoam.size() <<
endl;
402 if (tag !=
"$ENDNOD" && tag !=
"$EndNodes")
405 <<
"Did not find $ENDNOD tag on line "
415 Info<<
"Starting to read physical names at line " << inFile.
lineNumber()
425 Info<<
"Physical names:" << nNames <<
endl;
427 for (label i = 0; i < nNames; i++)
435 label nSpaces = lineStr.str().count(
' ');
441 Info<<
" " << regionI <<
'\t'
444 else if (nSpaces == 2)
451 Info<<
" " <<
"Line " << regionI <<
'\t'
454 else if (physType == 2)
456 Info<<
" " <<
"Surface " << regionI <<
'\t'
459 else if (physType == 3)
461 Info<<
" " <<
"Volume " << regionI <<
'\t'
477 if (tag !=
"$EndPhysicalNames")
480 <<
"Did not find $EndPhysicalNames tag on line "
488 const scalar versionFormat,
489 const bool keepOrientation,
549 for (label elemI = 0; elemI < nElems; elemI++)
555 label elmNumber, elmType, regPhys;
556 if (versionFormat >= 2)
558 lineStr >> elmNumber >> elmType;
567 for (label i = 1; i < nTags; i++)
576 label regElem, nNodes;
577 lineStr >> elmNumber >> elmType >> regPhys >> regElem >> nNodes;
581 if (elmType == MSHLINE)
584 lineStr >> meshPti >> meshPti;
586 else if (elmType == MSHTRI)
592 const auto regFnd = physToPatch.cfind(regPhys);
603 patchi = patchFaces.size();
605 patchFaces.setSize(patchi + 1);
606 patchToPhys.
setSize(patchi + 1);
608 Info<<
"Mapping region " << regPhys <<
" to Foam patch "
610 physToPatch.insert(regPhys, patchi);
611 patchToPhys[patchi] = regPhys;
617 else if (elmType == MSHQUAD)
620 >> quadPoints[0] >> quadPoints[1] >> quadPoints[2]
625 const auto regFnd = physToPatch.cfind(regPhys);
636 patchi = patchFaces.size();
638 patchFaces.setSize(patchi + 1);
639 patchToPhys.
setSize(patchi + 1);
641 Info<<
"Mapping region " << regPhys <<
" to Foam patch "
643 physToPatch.insert(regPhys, patchi);
644 patchToPhys[patchi] = regPhys;
648 patchFaces[patchi].
append(quadPoints);
650 else if (elmType == MSHTET)
671 else if (elmType == MSHPYR)
683 >> pyrPoints[0] >> pyrPoints[1] >> pyrPoints[2]
684 >> pyrPoints[3] >> pyrPoints[4];
688 cells[celli++].reset(pyr, pyrPoints);
692 else if (elmType == MSHPRISM)
704 >> prismPoints[0] >> prismPoints[1] >> prismPoints[2]
705 >> prismPoints[3] >> prismPoints[4] >> prismPoints[5];
709 cells[celli].reset(prism, prismPoints);
713 if (!keepOrientation && !correctOrientation(
points,
cell))
715 Info<<
"Inverting prism " << celli <<
endl;
717 prismPoints[0] =
cell[0];
718 prismPoints[1] =
cell[2];
719 prismPoints[2] =
cell[1];
720 prismPoints[3] =
cell[3];
721 prismPoints[4] =
cell[4];
722 prismPoints[5] =
cell[5];
724 cells[celli].reset(prism, prismPoints);
731 else if (elmType == MSHHEX)
743 >> hexPoints[0] >> hexPoints[1]
744 >> hexPoints[2] >> hexPoints[3]
745 >> hexPoints[4] >> hexPoints[5]
746 >> hexPoints[6] >> hexPoints[7];
754 if (!keepOrientation && !correctOrientation(
points,
cell))
756 Info<<
"Inverting hex " << celli <<
endl;
758 hexPoints[0] =
cell[4];
759 hexPoints[1] =
cell[5];
760 hexPoints[2] =
cell[6];
761 hexPoints[3] =
cell[7];
762 hexPoints[4] =
cell[0];
763 hexPoints[5] =
cell[1];
764 hexPoints[6] =
cell[2];
765 hexPoints[7] =
cell[3];
776 Info<<
"Unhandled element " << elmType <<
" at line "
786 if (tag !=
"$ENDELM" && tag !=
"$EndElements")
789 <<
"Did not find $ENDELM tag on line "
796 forAll(patchFaces, patchi)
798 patchFaces[patchi].shrink();
804 <<
" hex :" << nHex <<
endl
805 <<
" prism:" << nPrism <<
endl
806 <<
" pyr :" << nPyr <<
endl
807 <<
" tet :" << nTet <<
endl
810 if (
cells.size() == 0)
813 <<
"No cells read from file " << inFile.
name() <<
nl
814 <<
"Does your file specify any 3D elements (hex=" << MSHHEX
815 <<
", prism=" << MSHPRISM <<
", pyramid=" << MSHPYR
816 <<
", tet=" << MSHTET <<
")?" <<
nl
817 <<
"Perhaps you have not exported the 3D elements?"
822 <<
"Zone\tSize" <<
endl;
826 zoneCells[zonei].shrink();
828 const labelList& zCells = zoneCells[zonei];
832 Info<<
" " << zonei <<
'\t' << zCells.size() <<
endl;
840 const scalar versionFormat,
841 const bool keepOrientation,
875 label nEntities, nElems, minElemTag, maxElemTag;
876 lineStr >> nEntities >> nElems >> minElemTag >> maxElemTag;
897 for (label entityi = 0; entityi < nEntities; entityi++)
903 label entityDim, entityID, regPhys, elmType, nElemInBlock, elemID;
904 lineStr >> entityDim >> entityID >> elmType >> nElemInBlock;
907 regPhys = surfEntityToPhysSurface[entityID];
908 else if (entityDim == 3)
909 regPhys = volEntityToPhysVolume[entityID];
914 if (elmType == MSHLINE)
916 for (label entityElm = 0; entityElm < nElemInBlock; entityElm++)
921 lineStr >> elemID >> meshPti >> meshPti;
924 else if (elmType == MSHTRI)
926 for (label entityElm = 0; entityElm < nElemInBlock; entityElm++)
934 const auto regFnd = physToPatch.cfind(regPhys);
945 patchi = patchFaces.size();
947 patchFaces.setSize(patchi + 1);
948 patchToPhys.
setSize(patchi + 1);
950 Info<<
"Mapping region " << regPhys <<
" to Foam patch "
952 physToPatch.insert(regPhys, patchi);
953 patchToPhys[patchi] = regPhys;
960 else if (elmType == MSHQUAD)
962 for (label entityElm = 0; entityElm < nElemInBlock; entityElm++)
967 >> quadPoints[0] >> quadPoints[1] >> quadPoints[2]
972 const auto regFnd = physToPatch.cfind(regPhys);
983 patchi = patchFaces.size();
985 patchFaces.setSize(patchi + 1);
986 patchToPhys.
setSize(patchi + 1);
988 Info<<
"Mapping region " << regPhys <<
" to Foam patch "
990 physToPatch.insert(regPhys, patchi);
991 patchToPhys[patchi] = regPhys;
995 patchFaces[patchi].
append(quadPoints);
998 else if (elmType == MSHTET)
1000 nTet += nElemInBlock;
1002 for (label entityElm = 0; entityElm < nElemInBlock; entityElm++)
1025 else if (elmType == MSHPYR)
1027 nPyr += nElemInBlock;
1029 for (label entityElm = 0; entityElm < nElemInBlock; entityElm++)
1044 >> pyrPoints[0] >> pyrPoints[1] >> pyrPoints[2]
1045 >> pyrPoints[3] >> pyrPoints[4];
1049 cells[celli++].reset(pyr, pyrPoints);
1052 else if (elmType == MSHPRISM)
1054 nPrism += nElemInBlock;
1056 for (label entityElm = 0; entityElm < nElemInBlock; entityElm++)
1071 >> prismPoints[0] >> prismPoints[1] >> prismPoints[2]
1072 >> prismPoints[3] >> prismPoints[4] >> prismPoints[5];
1076 cells[celli].reset(prism, prismPoints);
1080 if (!keepOrientation && !correctOrientation(
points,
cell))
1082 Info<<
"Inverting prism " << celli <<
endl;
1084 prismPoints[0] =
cell[0];
1085 prismPoints[1] =
cell[2];
1086 prismPoints[2] =
cell[1];
1087 prismPoints[3] =
cell[3];
1088 prismPoints[4] =
cell[4];
1089 prismPoints[5] =
cell[5];
1091 cells[celli].reset(prism, prismPoints);
1097 else if (elmType == MSHHEX)
1099 nHex += nElemInBlock;
1101 for (label entityElm = 0; entityElm < nElemInBlock; entityElm++)
1116 >> hexPoints[0] >> hexPoints[1]
1117 >> hexPoints[2] >> hexPoints[3]
1118 >> hexPoints[4] >> hexPoints[5]
1119 >> hexPoints[6] >> hexPoints[7];
1123 cells[celli].reset(
hex, hexPoints);
1127 if (!keepOrientation && !correctOrientation(
points,
cell))
1129 Info<<
"Inverting hex " << celli <<
endl;
1131 hexPoints[0] =
cell[4];
1132 hexPoints[1] =
cell[5];
1133 hexPoints[2] =
cell[6];
1134 hexPoints[3] =
cell[7];
1135 hexPoints[4] =
cell[0];
1136 hexPoints[5] =
cell[1];
1137 hexPoints[6] =
cell[2];
1138 hexPoints[7] =
cell[3];
1140 cells[celli].reset(
hex, hexPoints);
1148 Info<<
"Unhandled element " << elmType <<
" at line "
1149 << inFile.
lineNumber() <<
"in/on physical region ID: "
1151 Info <<
"Perhaps you created a higher order mesh?" <<
endl;
1160 if (tag !=
"$ENDELM" && tag !=
"$EndElements")
1163 <<
"Did not find $ENDELM tag on line "
1170 forAll(patchFaces, patchi)
1172 patchFaces[patchi].shrink();
1178 <<
" hex : " << nHex <<
endl
1179 <<
" prism: " << nPrism <<
endl
1180 <<
" pyr : " << nPyr <<
endl
1181 <<
" tet : " << nTet <<
endl
1184 if (
cells.size() == 0)
1187 <<
"No cells read from file " << inFile.
name() <<
nl
1188 <<
"Does your file specify any 3D elements (hex=" << MSHHEX
1189 <<
", prism=" << MSHPRISM <<
", pyramid=" << MSHPYR
1190 <<
", tet=" << MSHTET <<
")?" <<
nl
1191 <<
"Perhaps you have not exported the 3D elements?"
1195 Info<<
"CellZones:" <<
nl
1196 <<
"Zone\tSize" <<
endl;
1200 zoneCells[zonei].shrink();
1202 const labelList& zCells = zoneCells[zonei];
1206 Info<<
" " << zonei <<
'\t' << zCells.size() <<
endl;
1219 label
nPoints, nCurves, nSurfaces, nVolumes;
1220 label entityID, physicalID, nPhysicalTags;
1226 lineStr >>
nPoints >> nCurves >> nSurfaces >> nVolumes;
1229 for (label i = 0; i <
nPoints; ++i)
1233 for (label i = 0; i < nCurves; ++i)
1237 for (label i = 0; i < nSurfaces; ++i)
1241 lineStr >> entityID;
1244 for (label j = 0; j < 6; ++j)
1248 lineStr >> nPhysicalTags;
1249 if (nPhysicalTags > 1)
1252 <<
"Cannot interpret multiple physical surfaces associated"
1253 <<
" with one surface on line number " << inFile.
lineNumber()
1257 lineStr >> physicalID;
1258 surfEntityToPhysSurface.insert(entityID, physicalID);
1262 for (label i = 0; i < nVolumes; ++i)
1266 lineStr >> entityID;
1269 for (label j = 0; j < 6; ++j)
1273 lineStr >> nPhysicalTags;
1274 if (nPhysicalTags > 1)
1277 <<
"Cannot interpret multiple physical volumes associated"
1278 <<
" with one volume on line number " << inFile.
lineNumber()
1282 lineStr >> physicalID;
1283 volEntityToPhysVolume.insert(entityID, physicalID);
1291 if (tag !=
"$EndEntities")
1294 <<
"Did not find $EndEntities tag on line "
1300 int main(
int argc,
char *argv[])
1304 "Convert a gmsh .msh file to OpenFOAM"
1312 "Retain raw orientation for prisms/hexs"
1327 const bool keepOrientation =
args.
found(
"keepOrientation");
1355 scalar versionFormat = 1;
1369 if (tag ==
"$MeshFormat")
1371 versionFormat = readMeshFormat(inFile);
1373 else if (tag ==
"$PhysicalNames")
1375 readPhysNames(inFile, physicalNames);
1377 else if (tag ==
"$Entities")
1380 readEntities(inFile,
1381 surfEntityToPhysSurface,
1382 volEntityToPhysVolume);
1384 else if (tag ==
"$NOD" || tag ==
"$Nodes")
1386 if (versionFormat < 4.0)
1387 readPointsLegacy(inFile,
points, mshToFoam);
1389 readPoints(inFile,
points, mshToFoam);
1391 else if (tag ==
"$ELM" || tag ==
"$Elements")
1393 if (versionFormat < 4.0)
1420 surfEntityToPhysSurface,
1421 volEntityToPhysVolume
1426 Info<<
"Skipping tag " << tag <<
" at line "
1430 if (!skipSection(inFile))
1435 }
while (inFile.
good());
1438 label nValidCellZones = 0;
1442 if (zoneCells[zonei].size())
1461 wordList boundaryPatchNames(boundaryFaces.size());
1463 forAll(boundaryPatchNames, patchi)
1465 boundaryPatchNames[patchi] =
1466 physicalNames.lookup
1468 patchToPhys[patchi],
1472 Info<<
"Patch " << patchi <<
" gets name "
1473 << boundaryPatchNames[patchi] <<
endl;
1477 wordList boundaryPatchTypes(boundaryFaces.size(), polyPatch::typeName);
1480 wordList boundaryPatchPhysicalTypes
1482 boundaryFaces.size(),
1501 boundaryPatchPhysicalTypes
1519 forAll(patchFaces, patchi)
1523 Info<<
"Finding faces of patch " << patchi <<
endl;
1530 label patchFacei = findFace(pp,
f);
1532 if (patchFacei != -1)
1534 label meshFacei = pp.
start() + patchFacei;
1536 repatcher.changePatchID(meshFacei, patchi);
1542 label meshFacei = findInternalFace(
mesh,
f);
1544 if (meshFacei != -1)
1546 zoneFaces[patchi].append(meshFacei);
1551 <<
"Could not match gmsh face " <<
f
1552 <<
" to any of the interior or exterior faces"
1553 <<
" that share the same 0th point" <<
endl;
1561 label nValidFaceZones = 0;
1563 Info<<
"FaceZones:" <<
nl
1564 <<
"Zone\tSize" <<
endl;
1568 zoneFaces[zonei].shrink();
1570 const labelList& zFaces = zoneFaces[zonei];
1576 Info<<
" " << zonei <<
'\t' << zFaces.size() <<
endl;
1586 repatcher.repatch();
1595 if (nValidCellZones > 0)
1599 nValidCellZones = 0;
1603 if (zoneCells[zonei].size())
1607 physicalNames.lookup
1614 Info<<
"Writing zone " << zonei <<
" to cellZone "
1615 << zoneName <<
" and cellSet"
1633 if (nValidFaceZones > 0)
1637 nValidFaceZones = 0;
1641 if (zoneFaces[zonei].size())
1645 physicalNames.lookup
1652 Info<<
"Writing zone " << zonei <<
" to faceZone "
1653 << zoneName <<
" and faceSet"
1672 if (cz.size() || fz.size())
1682 label newPatchi = 0;
1685 if (patchi != defaultPatchID)
1689 newPatchPtrList[newPatchi] =
patch.clone
1700 repatcher.changePatches(newPatchPtrList);