78void Foam::snappyLayerDriver::dumpDisplacement
80 const fileName& prefix,
81 const indirectPrimitivePatch& pp,
82 const vectorField& patchDisp,
83 const List<extrudeMode>& extrudeStatus
86 OBJstream dispStr(prefix +
"_disp.obj");
87 Info<<
"Writing all displacements to " << dispStr.name() <<
endl;
89 forAll(patchDisp, patchPointi)
91 const point& pt = pp.localPoints()[patchPointi];
92 dispStr.write(
linePointRef(pt, pt + patchDisp[patchPointi]));
96 OBJstream illStr(prefix +
"_illegal.obj");
97 Info<<
"Writing invalid displacements to " << illStr.name() <<
endl;
99 forAll(patchDisp, patchPointi)
101 if (extrudeStatus[patchPointi] !=
EXTRUDE)
103 const point& pt = pp.localPoints()[patchPointi];
104 illStr.write(
linePointRef(pt, pt + patchDisp[patchPointi]));
119 forAll(pp.localFaces(), facei)
121 const face&
f = pp.localFaces()[facei];
126 faceFld[facei] += pointFld[
f[fp]];
128 faceFld[facei] /=
f.
size();
137void Foam::snappyLayerDriver::checkManifold
140 pointSet& nonManifoldPoints
144 fp.checkPointManifold(
false, &nonManifoldPoints);
151 const labelList& eFaces = edgeFaces[edgei];
153 if (eFaces.size() > 2)
155 const edge&
e = fp.edges()[edgei];
157 nonManifoldPoints.insert(fp.meshPoints()[
e[0]]);
158 nonManifoldPoints.insert(fp.meshPoints()[
e[1]]);
164void Foam::snappyLayerDriver::checkMeshManifold()
const
166 const fvMesh&
mesh = meshRefiner_.mesh();
168 Info<<
nl <<
"Checking mesh manifoldness ..." <<
endl;
170 pointSet nonManifoldPoints
192 label nNonManif =
returnReduce(nonManifoldPoints.size(), sumOp<label>());
196 Info<<
"Outside of mesh is multiply connected across edges or"
198 <<
"This is not a fatal error but might cause some unexpected"
199 <<
" behaviour." <<
nl
214bool Foam::snappyLayerDriver::unmarkExtrusion
216 const label patchPointi,
219 List<extrudeMode>& extrudeStatus
222 if (extrudeStatus[patchPointi] == EXTRUDE)
224 extrudeStatus[patchPointi] = NOEXTRUDE;
225 patchNLayers[patchPointi] = 0;
226 patchDisp[patchPointi] =
Zero;
229 else if (extrudeStatus[patchPointi] == EXTRUDEREMOVE)
231 extrudeStatus[patchPointi] = NOEXTRUDE;
232 patchNLayers[patchPointi] = 0;
233 patchDisp[patchPointi] =
Zero;
242bool Foam::snappyLayerDriver::unmarkExtrusion
244 const face& localFace,
247 List<extrudeMode>& extrudeStatus
250 bool unextruded =
false;
272Foam::label Foam::snappyLayerDriver::constrainFp(
const label sz,
const label fp)
289void Foam::snappyLayerDriver::countCommonPoints
294 Map<label>& nCommonPoints
297 const faceList& localFaces = pp.localFaces();
300 const face&
f = localFaces[facei];
302 nCommonPoints.
clear();
306 label pointi =
f[fp];
311 label nbFacei =
pFaces[pFacei];
316 ++(nCommonPoints(nbFacei, 0));
323bool Foam::snappyLayerDriver::checkCommonOrder
333 const label nb = nbFace.find(curFace[fp]);
347 label fpPlus1 = curFace.fcIndex(fp);
348 label fpMin1 = curFace.rcIndex(fp);
351 label nbPlus1 = nbFace.fcIndex(nb);
352 label nbMin1 = nbFace.rcIndex(nb);
359 if (nbFace[nbPlus1] == curFace[fpPlus1])
364 else if (nbFace[nbPlus1] == curFace[fpMin1])
369 else if (nbFace[nbMin1] == curFace[fpMin1])
387 curFp = constrainFp(curFace.size(), curFp+curInc);
388 curNb = constrainFp(nbFace.size(), curNb+nbInc);
389 }
while (curFace[curFp] == nbFace[curNb]);
397 for (label commonI = 0; commonI < nCommon; commonI++)
399 curFp = constrainFp(curFace.size(), curFp+curInc);
400 curNb = constrainFp(nbFace.size(), curNb+nbInc);
402 if (curFace[curFp] != nbFace[curNb])
418void Foam::snappyLayerDriver::checkCommonOrder
422 const Map<label>& nCommonPoints,
425 List<extrudeMode>& extrudeStatus
430 const label nbFacei = iter.key();
431 const label nCommon = iter.val();
433 const face& curFace = pp[facei];
434 const face& nbFace = pp[nbFacei];
439 && nCommon != nbFace.size()
440 && nCommon != curFace.size()
443 bool stringOk = checkCommonOrder(nCommon, curFace, nbFace);
451 pp.localFaces()[facei],
458 pp.localFaces()[nbFacei],
469void Foam::snappyLayerDriver::handleNonStringConnected
474 List<extrudeMode>& extrudeStatus
482 List<extrudeMode> oldExtrudeStatus;
483 autoPtr<OBJstream> str;
486 oldExtrudeStatus = extrudeStatus;
491 meshRefiner_.mesh().time().path()
492 /
"nonStringConnected.obj"
495 Pout<<
"Dumping string edges to " << str().
name();
500 Map<label> nCommonPoints(128);
504 countCommonPoints(pp, facei, nCommonPoints);
527 forAll(extrudeStatus, pointi)
529 if (extrudeStatus[pointi] != oldExtrudeStatus[pointi])
533 meshRefiner_.mesh().points()[pp.meshPoints()[pointi]]
542void Foam::snappyLayerDriver::handleNonManifolds
549 List<extrudeMode>& extrudeStatus
552 const fvMesh&
mesh = meshRefiner_.mesh();
554 Info<<
nl <<
"Handling non-manifold points ..." <<
endl;
557 Info<<
nl <<
"Checking patch manifoldness ..." <<
endl;
559 pointSet nonManifoldPoints(
mesh,
"nonManifoldPoints", pp.
nPoints());
572 const labelList& eFaces = edgeFaces[edgei];
573 if (eFaces.size() > 2)
575 const edge&
e = pp.edges()[edgei];
576 nonManifoldPoints.insert(pp.meshPoints()[
e[0]]);
577 nonManifoldPoints.insert(pp.meshPoints()[
e[1]]);
583 forAll(edgeGlobalFaces, edgei)
585 if (edgeGlobalFaces[edgei].size() > 2)
590 const edge&
e = pp.edges()[edgei];
591 nonManifoldPoints.insert(pp.meshPoints()[
e[0]]);
592 nonManifoldPoints.insert(pp.meshPoints()[
e[1]]);
597 label nNonManif =
returnReduce(nonManifoldPoints.size(), sumOp<label>());
599 Info<<
"Outside of local patch is multiply connected across edges or"
600 <<
" points at " << nNonManif <<
" points." <<
endl;
607 nonManifoldPoints.sync(
mesh);
609 const labelList& meshPoints = pp.meshPoints();
611 forAll(meshPoints, patchPointi)
613 if (nonManifoldPoints.found(meshPoints[patchPointi]))
626 Info<<
"Set displacement to zero for all " << nNonManif
627 <<
" non-manifold points" <<
endl;
637 label nBaffleFaces = 0;
642 const labelList& fEdges = faceEdges[facei];
644 const labelList& globFaces0 = edgeGlobalFaces[fEdges[0]];
645 if (globFaces0.size() == 2)
647 const edge e0(globFaces0[0], globFaces0[1]);
648 bool isBaffle =
true;
649 for (label fp = 1; fp < fEdges.size(); fp++)
651 const labelList& globFaces = edgeGlobalFaces[fEdges[fp]];
654 (globFaces.size() != 2)
655 || (edge(globFaces[0], globFaces[1]) != e0)
665 bool unextrude = unmarkExtrusion
667 pp.localFaces()[facei],
685 reduce(nBaffleFaces, sumOp<label>());
689 Info<<
"Set displacement to zero for all points on " << nBaffleFaces
690 <<
" baffle faces" <<
endl;
697void Foam::snappyLayerDriver::handleFeatureAngle
701 const scalar minAngle,
704 List<extrudeMode>& extrudeStatus
707 const fvMesh&
mesh = meshRefiner_.mesh();
711 Info<<
nl <<
"Handling feature edges (angle < " << minAngle
714 if (minCos < 1-SMALL && minCos > -1+SMALL)
723 const labelList& eFaces = pp.edgeFaces()[edgei];
725 label meshEdgei = meshEdges[edgei];
731 edgeNormal[meshEdgei],
732 pp.faceNormals()[eFaces[i]]
745 autoPtr<OBJstream> str;
754 + meshRefiner_.timeName()
758 Info<<
"Writing feature edges to " << str().name() <<
endl;
767 const labelList& eFaces = pp.edgeFaces()[edgei];
769 label meshEdgei = meshEdges[edgei];
771 const vector&
n = edgeNormal[meshEdgei];
775 scalar
cos =
n & pp.faceNormals()[eFaces[0]];
779 const edge&
e = pp.edges()[edgei];
800 const point&
p0 = pp.localPoints()[
e[0]];
801 const point& p1 = pp.localPoints()[
e[1]];
808 Info<<
"Set displacement to zero for points on "
810 <<
" feature edges" <<
endl;
819void Foam::snappyLayerDriver::handleWarpedFaces
822 const scalar faceRatio,
824 const scalar edge0Len,
828 List<extrudeMode>& extrudeStatus
831 const fvMesh&
mesh = meshRefiner_.mesh();
834 Info<<
nl <<
"Handling cells with warped patch faces ..." <<
nl;
842 label nWarpedFaces = 0;
846 const face&
f = pp[i];
847 label faceI = pp.addressing()[i];
853 if (relativeSizes[patchI] &&
f.
size() > 3)
856 scalar edgeLen = edge0Len/(1<<ownLevel);
859 const point& fc = faceCentres[faceI];
860 const vector& fn = pp.faceNormals()[i];
867 vProj[fp] = (
n & fn);
871 scalar minVal =
min(vProj);
872 scalar maxVal =
max(vProj);
874 if ((maxVal - minVal) > faceRatio * edgeLen)
893 Info<<
"Set displacement to zero on "
895 <<
" warped faces since layer would be > " << faceRatio
896 <<
" of the size of the bounding box." <<
endl;
1003void Foam::snappyLayerDriver::setNumLayers
1009 List<extrudeMode>& extrudeStatus,
1013 const fvMesh&
mesh = meshRefiner_.mesh();
1015 Info<<
nl <<
"Handling points with inconsistent layer specification ..."
1025 label patchi = patchIDs[i];
1029 label wantedLayers = patchToNLayers[patchi];
1031 forAll(meshPoints, patchPointi)
1033 label ppPointi = pp.meshPointMap()[meshPoints[patchPointi]];
1035 maxLayers[ppPointi] =
max(wantedLayers, maxLayers[ppPointi]);
1036 minLayers[ppPointi] =
min(wantedLayers, minLayers[ppPointi]);
1067 <<
"Patchpoint:" << i <<
" coord:" << pp.localPoints()[i]
1068 <<
" maxLayers:" << maxLayers
1069 <<
" minLayers:" << minLayers
1072 else if (maxLayers[i] == minLayers[i])
1075 patchNLayers[i] = maxLayers[i];
1093 patchNLayers[i] = maxLayers[i];
1100 forAll(pp.localFaces(), facei)
1102 const face&
f = pp.localFaces()[facei];
1108 nCells =
max(nCells, patchNLayers[
f[fp]]);
1111 nAddedCells += nCells;
1113 reduce(nAddedCells, sumOp<label>());
1126Foam::snappyLayerDriver::makeLayerDisplacementField
1128 const pointMesh& pMesh,
1133 const pointBoundaryMesh& pointPatches = pMesh.boundary();
1137 pointPatches.size(),
1138 slipPointPatchVectorField::typeName
1140 wordList actualPatchTypes(patchFieldTypes.size());
1141 forAll(pointPatches, patchi)
1143 actualPatchTypes[patchi] = pointPatches[patchi].type();
1146 forAll(numLayers, patchi)
1150 if (numLayers[patchi] == 0)
1152 patchFieldTypes[patchi] =
1153 zeroFixedValuePointPatchVectorField::typeName;
1155 else if (numLayers[patchi] > 0)
1157 patchFieldTypes[patchi] = fixedValuePointPatchVectorField::typeName;
1161 forAll(pointPatches, patchi)
1163 if (isA<processorPointPatch>(pointPatches[patchi]))
1165 patchFieldTypes[patchi] = calculatedPointPatchVectorField::typeName;
1167 else if (isA<cyclicPointPatch>(pointPatches[patchi]))
1169 patchFieldTypes[patchi] = cyclicSlipPointPatchVectorField::typeName;
1174 const polyMesh&
mesh = pMesh();
1183 "pointDisplacement",
1197void Foam::snappyLayerDriver::growNoExtrusion
1202 List<extrudeMode>& extrudeStatus
1205 Info<<
nl <<
"Growing non-extrusion points by one layer ..." <<
endl;
1207 List<extrudeMode> grownExtrudeStatus(extrudeStatus);
1209 const faceList& localFaces = pp.localFaces();
1213 forAll(localFaces, facei)
1215 const face&
f = localFaces[facei];
1217 bool hasSqueeze =
false;
1220 if (extrudeStatus[
f[fp]] == NOEXTRUDE)
1234 extrudeStatus[
f[fp]] == EXTRUDE
1235 && grownExtrudeStatus[
f[fp]] != NOEXTRUDE
1238 grownExtrudeStatus[
f[fp]] = NOEXTRUDE;
1245 extrudeStatus.transfer(grownExtrudeStatus);
1254 status[i] = extrudeStatus[i];
1258 meshRefiner_.mesh(),
1266 extrudeStatus[i] = extrudeMode(status[i]);
1271 forAll(extrudeStatus, patchPointi)
1273 if (extrudeStatus[patchPointi] == NOEXTRUDE)
1275 patchDisp[patchPointi] =
Zero;
1276 patchNLayers[patchPointi] = 0;
1280 reduce(nGrown, sumOp<label>());
1282 Info<<
"Set displacement to zero for an additional " << nGrown
1283 <<
" points." <<
endl;
1287void Foam::snappyLayerDriver::determineSidePatches
1289 const globalIndex& globalFaces,
1305 fvMesh&
mesh = meshRefiner_.mesh();
1314 Map<label> nbrProcToPatch;
1315 Map<label> patchToNbrProc;
1336 Info<<
nl <<
"Adding in total " << nAdded/2 <<
" inter-processor patches to"
1337 <<
" handle extrusion of non-manifold processor boundaries."
1344 Map<label> wantedToAddedPatch;
1346 for (label patchi = nOldPatches; patchi <
nPatches; patchi++)
1348 label nbrProci = patchToNbrProc[patchi];
1354 dictionary patchDict;
1357 patchDict.add(
"neighbProcNo", nbrProci);
1358 patchDict.add(
"nFaces", 0);
1366 label procPatchi = meshRefiner_.appendPatch
1373 wantedToAddedPatch.insert(patchi, procPatchi);
1379 label patchi = edgePatchID[i];
1380 const auto fnd = wantedToAddedPatch.cfind(patchi);
1383 edgePatchID[i] = fnd.val();
1393void Foam::snappyLayerDriver::calculateLayerThickness
1397 const layerParameters& layerParams,
1400 const scalar edge0Len,
1407 const fvMesh&
mesh = meshRefiner_.mesh();
1415 scalarField firstLayerThickness(pp.nPoints(), GREAT);
1416 scalarField finalLayerThickness(pp.nPoints(), GREAT);
1420 minThickness.setSize(pp.nPoints());
1421 minThickness = GREAT;
1423 thickness.setSize(pp.nPoints());
1426 expansionRatio.setSize(pp.nPoints());
1427 expansionRatio = GREAT;
1429 for (
const label patchi : patchIDs)
1433 forAll(meshPoints, patchPointi)
1435 label ppPointi = pp.meshPointMap()[meshPoints[patchPointi]];
1437 firstLayerThickness[ppPointi] =
min
1439 firstLayerThickness[ppPointi],
1440 layerParams.firstLayerThickness()[patchi]
1442 finalLayerThickness[ppPointi] =
min
1444 finalLayerThickness[ppPointi],
1445 layerParams.finalLayerThickness()[patchi]
1447 totalThickness[ppPointi] =
min
1449 totalThickness[ppPointi],
1450 layerParams.thickness()[patchi]
1452 expRatio[ppPointi] =
min
1455 layerParams.expansionRatio()[patchi]
1457 minThickness[ppPointi] =
min
1459 minThickness[ppPointi],
1460 layerParams.minThickness()[patchi]
1469 firstLayerThickness,
1477 finalLayerThickness,
1517 label ownLevel = cellLevel[
mesh.
faceOwner()[pp.addressing()[i]]];
1519 const face&
f = pp.localFaces()[i];
1523 maxPointLevel[
f[fp]] =
max(maxPointLevel[
f[fp]], ownLevel);
1551 for (
const label patchi : patchIDs)
1555 layerParams.layerModels()[patchi];
1556 const bool relSize = layerParams.relativeSizes()[patchi];
1558 for (
const label meshPointi : meshPoints)
1560 const label ppPointi = pp.meshPointMap()[meshPointi];
1565 edgeLen[ppPointi] =
min
1568 edge0Len/(1<<maxPointLevel[ppPointi])
1570 spec[ppPointi] =
max(spec[ppPointi], patchSpec);
1571 isRelativePoint[meshPointi] =
1572 isRelativePoint[meshPointi]
1597 orEqOp<unsigned int>(),
1604 forAll(pp.meshPoints(), pointi)
1606 const label meshPointi = pp.meshPoints()[pointi];
1614 finalLayerThickness[pointi] *= edgeLen[pointi];
1615 if (isRelativePoint[meshPointi])
1617 totalThickness[pointi] *= edgeLen[pointi];
1618 minThickness[pointi] *= edgeLen[pointi];
1621 else if (isRelativePoint[meshPointi])
1623 firstLayerThickness[pointi] *= edgeLen[pointi];
1624 finalLayerThickness[pointi] *= edgeLen[pointi];
1625 totalThickness[pointi] *= edgeLen[pointi];
1626 minThickness[pointi] *= edgeLen[pointi];
1629 thickness[pointi] =
min
1635 patchNLayers[pointi],
1636 firstLayerThickness[pointi],
1637 finalLayerThickness[pointi],
1638 totalThickness[pointi],
1642 expansionRatio[pointi] =
min
1644 expansionRatio[pointi],
1645 layerParameters::layerExpansionRatio
1648 patchNLayers[pointi],
1649 firstLayerThickness[pointi],
1650 finalLayerThickness[pointi],
1651 totalThickness[pointi],
1689 label maxPatchNameLen = 0;
1692 label patchi = patchIDs[i];
1694 maxPatchNameLen =
max(maxPatchNameLen, label(patchName.size()));
1698 <<
setf(ios_base::left) <<
setw(maxPatchNameLen) <<
"patch"
1699 <<
setw(0) <<
" faces layers avg thickness[m]" <<
nl
1700 <<
setf(ios_base::left) <<
setw(maxPatchNameLen) <<
" "
1701 <<
setw(0) <<
" near-wall overall" <<
nl
1702 <<
setf(ios_base::left) <<
setw(maxPatchNameLen) <<
"-----"
1703 <<
setw(0) <<
" ----- ------ --------- -------" <<
endl;
1710 label patchi = patchIDs[i];
1714 layerParams.layerModels()[patchi];
1716 scalar sumThickness = 0;
1717 scalar sumNearWallThickness = 0;
1718 label nMasterPoints = 0;
1720 forAll(meshPoints, patchPointi)
1722 label meshPointi = meshPoints[patchPointi];
1723 if (isMasterPoint[meshPointi])
1725 label ppPointi = pp.meshPointMap()[meshPointi];
1727 sumThickness += thickness[ppPointi];
1728 sumNearWallThickness += layerParams.firstLayerThickness
1731 patchNLayers[ppPointi],
1732 firstLayerThickness[ppPointi],
1733 finalLayerThickness[ppPointi],
1734 thickness[ppPointi],
1735 expansionRatio[ppPointi]
1741 label totNPoints =
returnReduce(nMasterPoints, sumOp<label>());
1744 scalar avgThickness = 0;
1745 scalar avgNearWallThickness = 0;
1752 avgNearWallThickness =
1761 <<
" " <<
setw(6) << layerParams.numLayers()[patchi]
1762 <<
" " <<
setw(8) << avgNearWallThickness
1763 <<
" " <<
setw(8) << avgThickness
1772void Foam::snappyLayerDriver::syncPatchDisplacement
1778 List<extrudeMode>& extrudeStatus
1781 const fvMesh&
mesh = meshRefiner_.mesh();
1782 const labelList& meshPoints = pp.meshPoints();
1784 label nChangedTotal = 0;
1796 minMagSqrEqOp<vector>(),
1803 if (
mag(patchDisp[i]) < minThickness[i])
1821 labelList syncPatchNLayers(patchNLayers);
1834 forAll(syncPatchNLayers, i)
1836 if (syncPatchNLayers[i] != patchNLayers[i])
1865 forAll(syncPatchNLayers, i)
1867 if (syncPatchNLayers[i] != patchNLayers[i])
1884 nChangedTotal += nChanged;
1903void Foam::snappyLayerDriver::getPatchDisplacement
1912 List<extrudeMode>& extrudeStatus
1915 Info<<
nl <<
"Determining displacement for added points"
1916 <<
" according to pointNormal ..." <<
endl;
1918 const fvMesh&
mesh = meshRefiner_.mesh();
1921 const pointField& localPoints = pp.localPoints();
1933 patchDisp = thickness*pointNormals;
1936 label nNoVisNormal = 0;
1937 label nExtrudeRemove = 0;
2030 forAll(pointNormals, patchPointi)
2032 label meshPointi = pp.meshPoints()[patchPointi];
2034 if (extrudeStatus[patchPointi] == NOEXTRUDE)
2037 patchNLayers[patchPointi] = 0;
2038 patchDisp[patchPointi] =
Zero;
2043 const vector&
n = pointNormals[patchPointi];
2049 Pout<<
"No valid normal for point " << meshPointi
2050 <<
' ' << pp.points()[meshPointi]
2051 <<
"; setting displacement to "
2052 << patchDisp[patchPointi]
2056 extrudeStatus[patchPointi] = EXTRUDEREMOVE;
2063 forAll(extrudeStatus, patchPointi)
2065 if (extrudeStatus[patchPointi] == EXTRUDEREMOVE)
2070 const labelList& pEdges = pp.pointEdges()[patchPointi];
2074 label edgei = pEdges[i];
2076 label otherPointi = pp.edges()[edgei].otherVertex(patchPointi);
2078 if (extrudeStatus[otherPointi] != NOEXTRUDE)
2080 avg += localPoints[otherPointi] + patchDisp[otherPointi];
2089 Pout<<
"Displacement at illegal point "
2090 << localPoints[patchPointi]
2092 << (avg /
nPoints - localPoints[patchPointi])
2096 patchDisp[patchPointi] =
2098 - localPoints[patchPointi];
2111 <<
" points with point normal pointing through faces." <<
nl
2112 <<
"Reset displacement at "
2114 <<
" points to average of surrounding points." <<
endl;
2117 syncPatchDisplacement
2130bool Foam::snappyLayerDriver::sameEdgeNeighbour
2133 const label myGlobalFacei,
2134 const label nbrGlobFacei,
2138 const labelList& eFaces = globalEdgeFaces[edgei];
2139 if (eFaces.size() == 2)
2141 return edge(myGlobalFacei, nbrGlobFacei) == edge(eFaces[0], eFaces[1]);
2148void Foam::snappyLayerDriver::getVertexString
2154 const label myGlobFacei,
2155 const label nbrGlobFacei,
2159 const labelList& fEdges = pp.faceEdges()[facei];
2160 label fp = fEdges.
find(edgei);
2173 label prevFp = fEdges.rcIndex(startFp);
2193 label nextFp = fEdges.fcIndex(endFp);
2210 const face&
f = pp.localFaces()[facei];
2227Foam::label Foam::snappyLayerDriver::truncateDisplacement
2229 const globalIndex& globalFaces,
2233 const faceSet& illegalPatchFaces,
2236 List<extrudeMode>& extrudeStatus
2239 const fvMesh&
mesh = meshRefiner_.mesh();
2243 const Map<label>& meshPointMap = pp.meshPointMap();
2245 for (
const label facei : illegalPatchFaces)
2250 <<
"Faceset " << illegalPatchFaces.name()
2251 <<
" contains internal face " << facei <<
nl
2260 const auto fnd = meshPointMap.cfind(
f[fp]);
2263 const label patchPointi = fnd.val();
2265 if (extrudeStatus[patchPointi] != NOEXTRUDE)
2280 forAll(patchDisp, patchPointi)
2282 if (
mag(patchDisp[patchPointi]) < minThickness[patchPointi])
2298 else if (extrudeStatus[patchPointi] == NOEXTRUDE)
2301 patchDisp[patchPointi] =
Zero;
2302 patchNLayers[patchPointi] = 0;
2307 const faceList& localFaces = pp.localFaces();
2311 syncPatchDisplacement
2332 const face& localF = localFaces[i];
2337 extrudeMode prevMode = extrudeStatus[localF.prevLabel(0)];
2341 extrudeMode fpMode = extrudeStatus[localF[fp]];
2343 if (prevMode == NOEXTRUDE && fpMode != NOEXTRUDE)
2370 reduce(nPinched, sumOp<label>());
2372 Info<<
"truncateDisplacement : Unextruded " << nPinched
2373 <<
" faces due to non-consecutive vertices being extruded." <<
endl;
2396 label nButterFly = 0;
2398 DynamicList<label> stringedVerts;
2399 forAll(pp.edges(), edgei)
2401 const labelList& globFaces = edgeGlobalFaces[edgei];
2403 if (globFaces.size() == 2)
2405 label myFacei = pp.edgeFaces()[edgei][0];
2406 label myGlobalFacei = globalFaces.toGlobal
2408 pp.addressing()[myFacei]
2410 label nbrGlobalFacei =
2412 globFaces[0] != myGlobalFacei
2429 extrudeStatus[stringedVerts[0]] != NOEXTRUDE
2430 || extrudeStatus[stringedVerts.last()] != NOEXTRUDE
2435 for (label i = 1; i < stringedVerts.size()-1; i++)
2437 if (extrudeStatus[stringedVerts[i]] == NOEXTRUDE)
2468 reduce(nButterFly, sumOp<label>());
2470 Info<<
"truncateDisplacement : Unextruded " << nButterFly
2471 <<
" faces due to stringed edges with inconsistent extrusion."
2482 label nDiffering = 0;
2526 if (nPinched+nButterFly+nDiffering == 0)
2538void Foam::snappyLayerDriver::setupLayerInfoTruncation
2542 const List<extrudeMode>& extrudeStatus,
2543 const label nBufferCellsNoExtrude,
2548 Info<<
nl <<
"Setting up information for layer truncation ..." <<
endl;
2550 const fvMesh&
mesh = meshRefiner_.mesh();
2552 if (nBufferCellsNoExtrude < 0)
2554 Info<<
nl <<
"Performing no layer truncation."
2555 <<
" nBufferCellsNoExtrude set to less than 0 ..." <<
endl;
2558 forAll(pp.localFaces(), patchFacei)
2560 const face&
f = pp.localFaces()[patchFacei];
2564 const label nPointLayers = patchNLayers[
f[fp]];
2565 if (nPointLayers > 0)
2567 if (nPatchFaceLayers[patchFacei] == -1)
2569 nPatchFaceLayers[patchFacei] = nPointLayers;
2573 nPatchFaceLayers[patchFacei] =
min
2575 nPatchFaceLayers[patchFacei],
2582 nPatchPointLayers = patchNLayers;
2585 forAll(nPatchFaceLayers, patchFacei)
2587 if (nPatchFaceLayers[patchFacei] == -1)
2589 nPatchFaceLayers[patchFacei] = 0;
2598 forAll(pp.localFaces(), patchFacei)
2600 const face&
f = pp.localFaces()[patchFacei];
2605 bool noExtrude =
false;
2610 if (extrudeStatus[
f[fp]] == NOEXTRUDE)
2614 mLevel =
max(mLevel, patchNLayers[
f[fp]]);
2624 nPatchFaceLayers[patchFacei] = 1;
2625 maxLevel[patchFacei] = mLevel;
2629 maxLevel[patchFacei] = mLevel;
2641 label nLevels =
gMax(patchNLayers);
2644 for (label ilevel = 1; ilevel < nLevels; ilevel++)
2650 nBuffer = nBufferCellsNoExtrude - 1;
2654 nBuffer = nBufferCellsNoExtrude;
2657 for (label ibuffer = 0; ibuffer < nBuffer + 1; ibuffer++)
2659 labelList tempCounter(nPatchFaceLayers);
2661 boolList foundNeighbour(pp.nPoints(),
false);
2663 forAll(pp.meshPoints(), patchPointi)
2665 forAll(pointFaces[patchPointi], pointFacei)
2667 label facei = pointFaces[patchPointi][pointFacei];
2671 nPatchFaceLayers[facei] != -1
2672 && maxLevel[facei] > 0
2675 foundNeighbour[patchPointi] =
true;
2690 forAll(pp.meshPoints(), patchPointi)
2692 if (foundNeighbour[patchPointi])
2694 forAll(pointFaces[patchPointi], pointFacei)
2696 label facei = pointFaces[patchPointi][pointFacei];
2699 nPatchFaceLayers[facei] == -1
2700 && maxLevel[facei] > 0
2701 && ilevel < maxLevel[facei]
2704 tempCounter[facei] = ilevel;
2709 nPatchFaceLayers = tempCounter;
2713 forAll(pp.localFaces(), patchFacei)
2715 if (nPatchFaceLayers[patchFacei] == -1)
2717 nPatchFaceLayers[patchFacei] = maxLevel[patchFacei];
2721 forAll(pp.meshPoints(), patchPointi)
2723 if (extrudeStatus[patchPointi] != NOEXTRUDE)
2725 forAll(pointFaces[patchPointi], pointFacei)
2727 label face = pointFaces[patchPointi][pointFacei];
2728 nPatchPointLayers[patchPointi] =
max
2730 nPatchPointLayers[patchPointi],
2731 nPatchFaceLayers[face]
2737 nPatchPointLayers[patchPointi] = 0;
2753bool Foam::snappyLayerDriver::cellsUseFace
2755 const polyMesh&
mesh,
2762 const cell& cFaces =
mesh.
cells()[cellLabels[i]];
2766 if (faces.found(cFaces[cFacei]))
2780Foam::label Foam::snappyLayerDriver::checkAndUnmark
2782 const addPatchCellLayer& addLayer,
2783 const dictionary& meshQualityDict,
2784 const bool additionalReporting,
2785 const List<labelPair>& baffles,
2787 const fvMesh& newMesh,
2791 List<extrudeMode>& extrudeStatus
2795 Info<<
nl <<
"Checking mesh with layer ..." <<
endl;
2796 faceSet wrongFaces(newMesh,
"wrongFaces", newMesh.nFaces()/1000);
2809 <<
" (concave, zero area or negative cell pyramid volume)"
2823 addLayer.layerFaces()
2831 const label nReportMax = 10;
2832 DynamicField<point> disabledFaceCentres(nReportMax);
2834 forAll(addedCells, oldPatchFacei)
2838 const labelList& fCells = addedCells[oldPatchFacei];
2840 if (cellsUseFace(newMesh, fCells, wrongFaces))
2847 pp.localFaces()[oldPatchFacei],
2854 if (additionalReporting && (nChanged < nReportMax))
2856 disabledFaceCentres.
append
2858 pp.faceCentres()[oldPatchFacei]
2868 label nChangedTotal =
returnReduce(nChanged, sumOp<label>());
2870 if (additionalReporting)
2879 label nReportLocal = nChanged;
2880 if (nChangedTotal > nReportMax)
2895 Pout<<
"Checked mesh with layers. Disabled extrusion at " <<
endl;
2896 for (label i=0; i < nReportLocal; i++)
2898 Pout<<
" " << disabledFaceCentres[i] <<
endl;
2902 label nReportTotal =
returnReduce(nReportLocal, sumOp<label>());
2904 if (nReportTotal < nChangedTotal)
2906 Info<<
"Suppressed disabled extrusion message for other "
2907 << nChangedTotal - nReportTotal <<
" faces." <<
endl;
2911 return nChangedTotal;
2916Foam::label Foam::snappyLayerDriver::countExtrusion
2919 const List<extrudeMode>& extrudeStatus
2923 label nExtruded = 0;
2925 const faceList& localFaces = pp.localFaces();
2929 const face& localFace = localFaces[i];
2933 if (extrudeStatus[localFace[fp]] != NOEXTRUDE)
2948 const polyMesh&
mesh,
2950 const List<labelPair>& baffles
2959 Map<label> baffleSet(4*baffles.size());
2962 baffleSet.insert(baffles[bafflei][0], bafflei);
2963 baffleSet.insert(baffles[bafflei][1], bafflei);
2967 List<labelPair> newBaffles(baffles.size(),
labelPair(-1, -1));
2975 label oldFacei = newToOldFaces[facei];
2977 const auto faceFnd = baffleSet.find(oldFacei);
2978 if (faceFnd.found())
2980 label bafflei = faceFnd();
2986 else if (
p[1] == -1)
2993 <<
"Problem:" << facei <<
" at:"
2995 <<
" is on same baffle as " <<
p[0]
3008void Foam::snappyLayerDriver::getLayerCellsFaces
3010 const polyMesh&
mesh,
3011 const addPatchCellLayer& addLayer,
3021 faceRealThickness = 0;
3029 forAll(addedCells, oldPatchFacei)
3031 const labelList& added = addedCells[oldPatchFacei];
3033 const labelList& layer = layerFaces[oldPatchFacei];
3040 cellNLayers[added[i]] = layer.size()-1;
3045 forAll(layerFaces, oldPatchFacei)
3047 const labelList& layer = layerFaces[oldPatchFacei];
3048 const scalar realThickness = oldRealThickness[oldPatchFacei];
3054 for (label i = 1; i < layer.size(); i++)
3056 faceRealThickness[layer[i]] = realThickness;
3063void Foam::snappyLayerDriver::printLayerData
3077 label maxPatchNameLen = 0;
3080 label patchi = patchIDs[i];
3081 word patchName = pbm[patchi].name();
3082 maxPatchNameLen =
max(maxPatchNameLen, label(patchName.size()));
3086 <<
setf(ios_base::left) <<
setw(maxPatchNameLen) <<
"patch"
3087 <<
setw(0) <<
" faces layers overall thickness" <<
nl
3088 <<
setf(ios_base::left) <<
setw(maxPatchNameLen) <<
" "
3089 <<
setw(0) <<
" [m] [%]" <<
nl
3090 <<
setf(ios_base::left) <<
setw(maxPatchNameLen) <<
"-----"
3091 <<
setw(0) <<
" ----- ------ --- ---" <<
endl;
3096 label patchi = patchIDs[i];
3097 const polyPatch& pp = pbm[patchi];
3099 label sumSize = pp.size();
3102 const labelList& faceCells = pp.faceCells();
3103 label sumNLayers = 0;
3106 sumNLayers += cellNLayers[faceCells[i]];
3119 scalar sumRealThickness =
sum(patchReal);
3120 scalar sumFraction = 0;
3123 if (patchWanted[i] > VSMALL)
3125 sumFraction += (patchReal[i]/patchWanted[i]);
3130 reduce(sumSize, sumOp<label>());
3131 reduce(sumNLayers, sumOp<label>());
3132 reduce(sumRealThickness, sumOp<scalar>());
3133 reduce(sumFraction, sumOp<scalar>());
3136 scalar avgLayers = 0;
3138 scalar avgFraction = 0;
3141 avgLayers = scalar(sumNLayers)/sumSize;
3142 avgReal = sumRealThickness/sumSize;
3143 avgFraction = sumFraction/sumSize;
3148 <<
" " <<
setw(8) << sumSize
3149 <<
" " <<
setw(8) << avgLayers
3150 <<
" " <<
setw(8) << avgReal
3151 <<
" " <<
setw(8) << 100*avgFraction
3158bool Foam::snappyLayerDriver::writeLayerSets
3168 forAll(cellNLayers, celli)
3170 if (cellNLayers[celli] > 0)
3175 cellSet addedCellSet(
mesh,
"addedCells", nAdded);
3176 forAll(cellNLayers, celli)
3178 if (cellNLayers[celli] > 0)
3180 addedCellSet.insert(celli);
3183 addedCellSet.instance() = meshRefiner_.timeName();
3186 <<
" added cells to cellSet "
3187 << addedCellSet.name() <<
endl;
3188 bool ok = addedCellSet.
write();
3189 allOk = allOk && ok;
3195 if (faceRealThickness[facei] > 0)
3201 faceSet layerFacesSet(
mesh,
"layerFaces", nAdded);
3204 if (faceRealThickness[facei] > 0)
3206 layerFacesSet.insert(facei);
3209 layerFacesSet.instance() = meshRefiner_.timeName();
3212 <<
" faces inside added layer to faceSet "
3213 << layerFacesSet.name() <<
endl;
3214 bool ok = layerFacesSet.
write();
3215 allOk = allOk && ok;
3221bool Foam::snappyLayerDriver::writeLayerData
3234 bool ok = writeLayerSets(
mesh, cellNLayers, faceRealThickness);
3235 allOk = allOk && ok;
3256 fixedValueFvPatchScalarField::typeName
3260 fld[celli] = cellNLayers[celli];
3267 label patchi = patchIDs[i];
3268 const polyPatch& pp = pbm[patchi];
3269 const labelList& faceCells = pp.faceCells();
3273 pfld[i] = cellNLayers[faceCells[i]];
3275 fldBf[patchi] == pfld;
3279 bool ok =
fld.write();
3280 allOk = allOk && ok;
3296 fixedValueFvPatchScalarField::typeName
3302 label patchi = patchIDs[i];
3303 fldBf[patchi] == pbm[patchi].patchSlice(faceRealThickness);
3307 bool ok =
fld.write();
3308 allOk = allOk && ok;
3315 "thicknessFraction",
3324 fixedValueFvPatchScalarField::typeName
3330 label patchi = patchIDs[i];
3345 if (patchWanted[i] > VSMALL)
3347 pfld[i] = patchReal[i]/patchWanted[i];
3351 fldBf[patchi] == pfld;
3354 <<
" : overall layer thickness (fraction"
3355 <<
" of desired thickness)" <<
endl;
3356 bool ok =
fld.write();
3357 allOk = allOk && ok;
3366void Foam::snappyLayerDriver::dupFaceZonePoints
3370 List<labelPair> baffles,
3375 fvMesh&
mesh = meshRefiner_.mesh();
3384 autoPtr<indirectPrimitivePatch> pp
3398 label nIdealTotAddedCells = 0;
3399 List<extrudeMode> extrudeStatus(pp().
nPoints(), EXTRUDE);
3414 syncPatchDisplacement
3428 forAll(extrudeStatus, patchPointi)
3430 label pointi = pp().meshPoints()[patchPointi];
3431 minPatchState[pointi] = extrudeStatus[patchPointi];
3449 forAll(minPatchState, pointi)
3451 label state = minPatchState[pointi];
3452 if (state == EXTRUDE || state == EXTRUDEREMOVE)
3457 candidatePoints.setSize(
n);
3459 forAll(minPatchState, pointi)
3461 label state = minPatchState[pointi];
3462 if (state == EXTRUDE || state == EXTRUDEREMOVE)
3464 candidatePoints[
n++] = pointi;
3482 bool hasInfo = meshRefiner_.getFaceZoneInfo
3489 if (hasInfo && !layerIDs.found(mpi) && !layerIDs.found(spi))
3491 nonDupZones.append(zonei);
3503 const localPointRegion regionSide(
mesh, nonDupBaffles, candidatePoints);
3505 autoPtr<mapPolyMesh> map = meshRefiner_.dupNonManifoldPoints
3515 const labelList& pointMap = map().pointMap();
3516 const labelList& reversePointMap = map().reversePointMap();
3520 label oldPointi = pointMap[pointi];
3521 label newMasterPointi = reversePointMap[oldPointi];
3523 if (newMasterPointi != pointi)
3526 pointToMaster[pointi] = newMasterPointi;
3527 pointToMaster[newMasterPointi] = newMasterPointi;
3533 const labelList& reverseFaceMap = map().reverseFaceMap();
3536 label f0 = reverseFaceMap[baffles[i].first()];
3537 label f1 = reverseFaceMap[baffles[i].second()];
3546 Info<<
"Writing point-duplicate mesh to time "
3547 << meshRefiner_.timeName() <<
endl;
3563 /
"duplicatePoints_"
3564 + meshRefiner_.timeName()
3567 Info<<
"Writing point-duplicates to " << str.name() <<
endl;
3571 label newMasteri = reversePointMap[pointMap[pointi]];
3573 if (newMasteri != pointi)
3583void Foam::snappyLayerDriver::mergeFaceZonePoints
3594 fvMesh&
mesh = meshRefiner_.mesh();
3597 label nPointPairs = 0;
3598 forAll(pointToMaster, pointi)
3600 label otherPointi = pointToMaster[pointi];
3601 if (otherPointi != -1)
3606 reduce(nPointPairs, sumOp<label>());
3607 if (nPointPairs > 0)
3610 Info<<
"Merging " << nPointPairs <<
" duplicated points ..." <<
endl;
3618 + meshRefiner_.timeName()
3621 Info<<
"Points to be merged to " << str.name() <<
endl;
3622 forAll(pointToMaster, pointi)
3624 label otherPointi = pointToMaster[pointi];
3625 if (otherPointi != -1)
3635 autoPtr<mapPolyMesh> map = meshRefiner_.mergePoints(pointToMaster);
3640 const labelList& reverseFaceMap = map().reverseFaceMap();
3644 Info<<
"Merged points in = "
3652 Info<<
"Converting baffles back into zoned faces ..."
3655 autoPtr<mapPolyMesh> map = meshRefiner_.mergeZoneBaffles
3670 forAll(newFaceRealThickness, facei)
3672 label oldFacei =
faceMap[facei];
3675 scalar& realThick = newFaceRealThickness[facei];
3676 realThick =
max(realThick, faceRealThickness[oldFacei]);
3677 scalar& wanted = newFaceWantedThickness[facei];
3678 wanted =
max(wanted, faceWantedThickness[oldFacei]);
3681 faceRealThickness.transfer(newFaceRealThickness);
3682 faceWantedThickness.transfer(newFaceWantedThickness);
3685 Info<<
"Converted baffles in = "
3686 << meshRefiner_.mesh().time().cpuTimeIncrement()
3692Foam::label Foam::snappyLayerDriver::setPointNumLayers
3694 const layerParameters& layerParams,
3703 List<extrudeMode>& extrudeStatus
3706 fvMesh&
mesh = meshRefiner_.mesh();
3708 patchDisp.setSize(pp.nPoints());
3709 patchDisp =
vector(GREAT, GREAT, GREAT);
3713 patchNLayers.setSize(pp.nPoints());
3714 patchNLayers =
Zero;
3717 label nIdealTotAddedCells = 0;
3720 extrudeStatus.setSize(pp.nPoints());
3721 extrudeStatus = EXTRUDE;
3745 handleNonStringConnected
3775 layerParams.featureAngle(),
3786 if (!layerParams.relativeSizes().found(
false))
3789 const scalar edge0Len =
3790 meshRefiner_.meshCutter().level0EdgeLength();
3791 const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
3796 layerParams.maxFaceThicknessRatio(),
3797 layerParams.relativeSizes(),
3822 for (label i = 0; i < layerParams.nGrow(); i++)
3832 return nIdealTotAddedCells;
3837Foam::snappyLayerDriver::makeMeshMover
3839 const layerParameters& layerParams,
3840 const dictionary& motionDict,
3848 fvMesh&
mesh = meshRefiner_.mesh();
3852 dictionary combinedDict(layerParams.dict());
3854 combinedDict.merge(motionDict);
3856 combinedDict.add(
"minThicknessName", minThickness.name());
3858 const List<labelPair> internalBaffles
3870 autoPtr<Foam::externalDisplacementMeshMover> medialAxisMoverPtr
3874 layerParams.meshShrinker(),
3887 if (errorMsg.size() || IOerrorMsg.size())
3896 <<
"Missing/incorrect required dictionary entries:"
3898 << IOerrorMsg.c_str() <<
nl
3899 << errorMsg.c_str() <<
nl <<
nl
3900 <<
"Exiting dry-run" <<
nl <<
endl;
3904 Perr<<
"\nFOAM parallel run exiting\n" <<
endl;
3914 return medialAxisMoverPtr;
3921 const label nLayerIter,
3925 const label nRelaxedIter,
3926 const label nAllowableErrors,
3932 const label nIdealTotAddedCells,
3973 makeLayerDisplacementField
3997 for (label iteration = 0; iteration < nLayerIter; iteration++)
4000 <<
"Layer addition iteration " << iteration <<
nl
4001 <<
"--------------------------" <<
endl;
4007 iteration < nRelaxedIter
4009 : motionDict.
subDict(
"relaxed")
4012 if (iteration >= nRelaxedIter)
4014 Info<<
"Switched to relaxed meshQuality constraints." <<
endl;
4023 syncPatchDisplacement
4033 getPatchDisplacement
4054 patchDisp = -patchDisp;
4073 combinedDict.
merge(motionDict);
4075 combinedDict.
merge(meshQualityDict);
4077 combinedDict.
add(
"minThicknessName", minThickness.
name());
4080 medialAxisMoverPtr().move
4100 truncateDisplacement
4125 Info<<
"Writing shrunk mesh to time "
4126 << meshRefiner_.timeName() <<
endl;
4155 labelList nPatchFaceLayers(pp.size(), -1);
4156 setupLayerInfoTruncation
4170 forAll(nPatchPointLayers, i)
4174 nPatchPointLayers[i],
4177 finalDisp[i] = ratio*patchDisp[i];
4181 const scalarField invExpansionRatio(1.0/expansionRatio);
4216 savedMeshMod = meshMod;
4238 fvMesh& newMesh = *newMeshPtr;
4242 addProfiling(grow,
"snappyHexMesh::layers::updateMesh");
4275 avgPointData(pp,
mag(patchDisp))(),
4283 label nAddedCells = 0;
4284 forAll(cellNLayers, celli)
4286 if (cellNLayers[celli] > 0)
4295 Info<<
"Writing layer mesh to time " << meshRefiner_.timeName()
4298 writeLayerSets(newMesh, cellNLayers, faceRealThickness);
4317 facei < newMesh.
nFaces();
4321 label newMeshFacei = map.
faceMap()[facei];
4322 if (newMeshFacei != -1)
4324 meshToNewMesh[newMeshFacei] = facei;
4335 meshToNewMesh[
p[0]],
4338 if (newMeshBaffle[0] != -1 && newMeshBaffle[1] != -1)
4340 newMeshBaffles[newi++] = newMeshBaffle;
4354 <<
" baffles across faceZones of type internal" <<
nl
4358 label nTotChanged = checkAndUnmark
4372 label nTotExtruded = countExtrusion(pp, extrudeStatus);
4376 Info<<
"Extruding " << nTotExtruded
4377 <<
" out of " << nTotFaces
4378 <<
" faces (" << 100.0*nTotExtruded/nTotFaces <<
"%)."
4379 <<
" Removed extrusion at " << nTotChanged <<
" faces."
4381 <<
"Added " << nTotAddedCells <<
" out of "
4382 << nIdealTotAddedCells
4383 <<
" cells (" << 100.0*nTotAddedCells/nIdealTotAddedCells
4386 if (nTotChanged == 0)
4394 medialAxisMoverPtr().movePoints(
mesh.
points());
4400 for (label i = 0; i < layerParams.
nGrow(); i++)
4416void Foam::snappyLayerDriver::mapFaceZonePoints
4444 oldFaceToBaffle.insert(baffle[0], i);
4445 oldFaceToBaffle.insert(baffle[1], i);
4457 const label oldFacei = map.
faceMap()[facei];
4458 const auto iter = oldFaceToBaffle.
find(oldFacei);
4459 if (oldFacei != -1 && iter.found())
4461 const label bafflei = iter();
4462 auto& newBaffle = newBaffles[bafflei];
4463 if (newBaffle[0] == -1)
4465 newBaffle[0] = facei;
4467 else if (newBaffle[1] == -1)
4469 newBaffle[1] = facei;
4475 <<
" already maps to baffle faces:"
4478 <<
" and " << newBaffle[1]
4486 label pointi =
f[fp];
4487 label oldPointi = map.
pointMap()[pointi];
4489 if (pointToMaster[oldPointi] != -1)
4491 candidates.
append(pointi);
4503 const labelPair& newBaffle = newBaffles[i];
4504 if (newBaffle[0] != -1 && newBaffle[1] != -1)
4506 newBaffles[
n++] = newBaffle;
4524 UIndirectList<point>(
mesh.
points(), candidates),
4525 meshRefiner_.mergeDistance(),
4541 const labelList& oldPoints = newToOld[newi];
4542 if (oldPoints.size() > 1)
4548 label masteri =
min(meshPoints);
4551 pointToMaster[meshPoints[i]] = masteri;
4559void Foam::snappyLayerDriver::updatePatch
4562 const mapPolyMesh& map,
4563 autoPtr<indirectPrimitivePatch>& pp,
4569 fvMesh&
mesh = meshRefiner_.mesh();
4571 autoPtr<indirectPrimitivePatch> newPp
4581 newToOldPatchPoints.setSize(newPp().
nPoints());
4582 newToOldPatchPoints = -1;
4584 const Map<label>& baseMap = pp().meshPointMap();
4585 const labelList& newMeshPoints = newPp().meshPoints();
4589 const label newMeshPointi = newMeshPoints[
newPointi];
4590 const label oldMeshPointi =
4591 map.pointMap()[newMeshPointi];
4592 const auto iter = baseMap.find(oldMeshPointi);
4595 newToOldPatchPoints[
newPointi] = iter();
4604 pp = std::move(newPp);
4607 (void)pp().meshPoints();
4608 (void)pp().meshPointMap();
4622 meshRefiner_(meshRefiner),
4623 globalToMasterPatch_(globalToMasterPatch),
4624 globalToSlavePatch_(globalToSlavePatch),
4646 <<
"Merging all faces of a cell" <<
nl
4647 <<
"---------------------------" <<
nl
4648 <<
" - which are on the same patch" <<
nl
4649 <<
" - which make an angle < " << planarAngle
4651 <<
" (cos:" << minCos <<
')' <<
nl
4652 <<
" - as long as the resulting face doesn't become concave"
4655 <<
" (0=straight, 180=fully concave)" <<
nl
4666 duplicateFace[cpl[0]] = cpl[1];
4667 duplicateFace[cpl[1]] = cpl[0];
4670 label nChanged = meshRefiner_.mergePatchFacesUndo
4674 meshRefiner_.meshedPatches(),
4680 nChanged += meshRefiner_.mergeEdgesUndo(minCos, motionDict);
4689 const label nAllowableErrors,
4697 const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
4705 internalOrBaffleFaceZones = meshRefiner_.getZones(fzTypes);
4712 meshRefiner_.getZones
4727 meshRefiner_.createZoneBaffles
4737 Info<<
"Writing baffled mesh to time "
4738 << meshRefiner_.timeName() <<
endl;
4759 const labelList& reverseFaceMap = map->reverseFaceMap();
4762 label f0 = reverseFaceMap[baffles[i].
first()];
4763 label f1 = reverseFaceMap[baffles[i].second()];
4782 bool hasInfo = meshRefiner_.getFaceZoneInfo
4792 if (layerIDs.
found(mpi) && !layerIDs.
found(spi))
4796 <<
" adding layers to master patch " << pbm[mpi].
name()
4797 <<
" only. Freezing points on slave patch "
4801 else if (!layerIDs.
found(mpi) && layerIDs.
found(spi))
4805 <<
" adding layers to slave patch " << pbm[spi].
name()
4806 <<
" only. Freezing points on master patch "
4846 (void)pp().meshPoints();
4847 (void)pp().meshPointMap();
4881 const label nIdealTotAddedCells = setPointNumLayers
4903 meshRefiner_.timeName(),
4911 calculateLayerThickness
4916 meshRefiner_.meshCutter().cellLevel(),
4934 avgPointData(pp(), baseThickness);
4942 (basePatchNLayers+layerParams.
nOuterIter()-1)
4950 for (label layeri = 0; layeri < layerParams.
nOuterIter(); layeri++)
4962 const label nToAdd =
gSum(deltaNLayers);
4965 Info<<
"Outer iteration : " << layeri
4966 <<
" to add in current iteration : " << nToAdd <<
endl;
4981 determineSidePatches
5018 forAll(baseThickness, pointi)
5022 basePatchNLayers[pointi],
5023 baseThickness[pointi],
5024 baseExpansionRatio[pointi],
5025 basePatchNLayers[pointi]
5026 -nAddedLayers[pointi]
5027 -patchNLayers[pointi],
5028 patchNLayers[pointi]
5056 nIdealTotAddedCells,
5082 sliceFaceRealThickness
5087 const label nTotalAdded =
gSum(patchNLayers);
5090 Info<<
"Outer iteration : " << layeri
5091 <<
" added in current iteration : " << nTotalAdded
5092 <<
" out of : " <<
gSum(deltaNLayers) <<
endl;
5094 if (nTotalAdded == 0)
5101 forAll(patchNLayers, pointi)
5103 nAddedLayers[pointi] += patchNLayers[pointi];
5105 if (patchNLayers[pointi] == 0)
5116 basePatchNLayers[pointi] = nAddedLayers[pointi];
5117 deltaNLayers[pointi] = 0;
5125 deltaNLayers[pointi] =
max
5130 deltaNLayers[pointi],
5131 basePatchNLayers[pointi] - nAddedLayers[pointi]
5168 meshRefiner_.updateMesh(map,
labelList(0));
5180 cellNLayers[i] += sliceCellNLayers[i];
5190 faceRealThickness += sliceFaceRealThickness;
5206 faceWantedThickness,
5215 Info<<
"Writing mesh with layers but disconnected to time "
5216 << meshRefiner_.timeName() <<
endl;
5231 mapFaceZonePoints(map, baffles, pointToMaster);
5235 updatePatch(patchIDs, map, pp, newToOldPatchPoints);
5251 newToOldPatchPoints,
5257 newToOldPatchPoints,
5263 newToOldPatchPoints,
5264 extrudeMode::NOEXTRUDE,
5269 newToOldPatchPoints,
5275 newToOldPatchPoints,
5281 newToOldPatchPoints,
5287 newToOldPatchPoints,
5293 newToOldPatchPoints,
5317 <<
"Doing final balancing" <<
nl
5318 <<
"---------------------" <<
nl
5337 map().distributeCellData(cellNLayers);
5338 map().distributeFaceData(faceWantedThickness);
5339 map().distributeFaceData(faceRealThickness);
5353 faceWantedThickness,
5366 const bool preBalance,
5375 <<
"Shrinking and layer addition phase" <<
nl
5376 <<
"----------------------------------" <<
nl
5380 Info<<
"Using mesh parameters " << motionDict <<
nl <<
endl;
5389 mergePatchFacesUndo(layerParams, motionDict, mergeType);
5398 label nFacesWithLayers = 0;
5399 forAll(numLayers, patchi)
5401 if (numLayers[patchi] > 0)
5413 <<
"Ignoring layers on coupled patch " << pp.
name()
5425 meshRefiner_.getFaceZoneInfo(fZones[zonei].
name(), mpi, spi, fzType);
5427 if (numLayers[mpi] > 0)
5429 nFacesWithLayers += fZones[zonei].
size();
5431 if (numLayers[spi] > 0)
5433 nFacesWithLayers += fZones[zonei].
size();
5442 Info<<
nl <<
"No layers to generate ..." <<
endl;
5447 checkMeshManifold();
5450 Info<<
"Checking initial mesh ..." <<
endl;
5459 Info<<
"Detected " << nInitErrors <<
" illegal faces"
5460 <<
" (concave, zero area or negative cell pyramid volume)"
5464 bool faceZoneOnCoupledFace =
false;
5478 const faceZone& fZone = fZones[zonei];
5483 meshRefiner_.getFaceZoneInfo(fzName, mpi, spi, fzType);
5485 if (numLayers[mpi] > 0 || numLayers[spi])
5487 isExtrudedZoneFace.
set(fZone);
5504 if (intOrCoupled[facei] && isExtrudedZoneFace.
test(facei))
5506 faceZoneOnCoupledFace =
true;
5521 <<
"Doing initial balancing" <<
nl
5522 <<
"-----------------------" <<
nl
5526 forAll(numLayers, patchi)
5528 if (numLayers[patchi] > 0)
5533 cellWeights[pp.
faceCells()[i]] += numLayers[patchi];
5542 const faceZone& fZone = fZones[zonei];
5547 meshRefiner_.getFaceZoneInfo(fzName, mpi, spi, fzType);
5549 if (numLayers[mpi] > 0)
5558 cellWeights[
cellIDs[i]] += numLayers[mpi];
5562 if (numLayers[spi] > 0)
5569 cellWeights[
cellIDs[i]] += numLayers[mpi];
Istream and Ostream manipulators taking arguments.
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
DynamicList< T, SizeMin > & shrink()
Shrink the allocated space to the number of elements used.
void append(const T &val)
Copy append an element to the end of this list.
SubField< scalar > subField
Declare type of subField.
void setSize(const label n)
Dummy function, to make FixedList consistent with List.
GeometricBoundaryField< scalar, fvPatchField, volMesh > Boundary
Type of boundary fields.
bool found(const Key &key) const
Return true if hashed entry is found in table.
label size() const noexcept
The number of elements in table.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
const word & name() const noexcept
Return the object name.
const fileName & instance() const noexcept
Read access to instance path component.
writeOption writeOpt() const noexcept
The write option.
void exit()
Job end with "exit" termination.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
void transfer(List< T > &list)
void setSize(const label n)
Alias for resize()
void append(const T &val)
Append an element at the end of the list.
void clear()
Clear the list, i.e. set size to zero.
A HashTable to objects of type <T> with a label key.
virtual const fileName & name() const
Get the name of the stream.
virtual int precision() const
Get precision of output field.
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
A list of faces which address into the list of points.
label nPoints() const
Number of points supporting patch faces.
const Field< point_type > & localPoints() const
Return pointField of points in patch.
virtual void movePoints(const Field< point_type > &)
Correct patch after moving points.
label nProcs() const noexcept
Number of ranks associated with PstreamBuffers.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
fileName path() const
Return path.
static word timeName(const scalar t, const int precision=precision_)
A List with indirect addressing. Like IndirectList but does not store addressing.
T & first()
Return the first element of the list.
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
void size(const label n)
Older name for setAddressableSize.
label fcIndex(const label i) const noexcept
static bool & parRun() noexcept
Test if this a parallel run.
label size() const noexcept
The number of elements in the list.
Adds layers of cells to outside of polyPatch. Can optionally create stand-alone extruded mesh (addToM...
void setRefinement(const globalIndex &globalFaces, const labelListList &globalEdgeFaces, const scalarField &expansionRatio, const indirectPrimitivePatch &pp, const bitSet &flip, const labelList &sidePatchID, const labelList &sideZoneID, const boolList &sideFlip, const labelList &inflateFaceID, const labelList &exposedPatchID, const labelList &nFaceLayers, const labelList &nPointLayers, const vectorField &firstLayerDisp, polyTopoChange &meshMod)
Play commands into polyTopoChange to create layers on top.
void updateMesh(const mapPolyMesh &, const labelList &faceMap, const labelList &pointMap)
Update any locally stored mesh information. Gets additional.
static void calcExtrudeInfo(const bool zoneFromAnyFace, const polyMesh &, const globalIndex &globalFaces, const labelListList &globalEdgeFaces, const indirectPrimitivePatch &pp, labelList &edgePatchID, label &nPatches, Map< label > &nbrProcToPatch, Map< label > &patchToNbrProc, labelList &edgeZoneID, boolList &edgeFlip, labelList &inflateFaceID)
Determine extrude information per patch edge:
labelListList addedCells() const
Added cells given current mesh & layerfaces.
static labelListList globalEdgeFaces(const polyMesh &, const globalIndex &globalFaces, const indirectPrimitivePatch &pp)
Per patch edge the pp faces (in global indices) using it.
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
void set(const bitSet &bitset)
Set specified bits from another bitset.
bool test(const label pos) const
Test value at specified position, never auto-vivify entries.
double cpuTimeIncrement() const
Return CPU time (in seconds) since last call to cpuTimeIncrement()
virtual const vectorField & pointNormals() const
Return point unit normals.
Abstract base class for domain decomposition.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
bool merge(const dictionary &dict)
Merge entries from the given dictionary.
entry * add(entry *entryPtr, bool mergeEntry=false)
Add a new entry.
string message() const
The accumulated error message.
A subset of mesh faces organised as a primitive patch.
const labelList & masterCells() const
const labelList & slaveCells() const
Return labels of slave cells.
Sends/receives parts of mesh+fvfields to neighbouring processors. Used in load balancing.
Mesh data needed to do the Finite Volume discretisation.
virtual bool write(const bool valid=true) const
Write mesh using IO settings from time.
const Time & time() const
Return the top-level database.
const word & name() const
Return reference to name.
virtual void updateMesh(const mapPolyMesh &mpm)
Update mesh corresponding to the given map.
virtual void movePoints(const pointField &)
Move points, returns volumes swept by faces in motion.
void clearOut()
Clear all geometry and addressing.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
void reset(const label localSize, const label comm=UPstream::worldComm, const bool parallel=UPstream::parRun())
void checkMesh() const
Debug: Check coupled mesh for correctness.
Simple container to keep together layer specific information.
scalar mergePatchFacesAngle() const
label nRelaxedIter() const
Number of iterations after which relaxed motion rules.
bool additionalReporting() const
Any additional reporting requested?
thicknessModelType
Enumeration defining the layer specification:
@ FIRST_AND_RELATIVE_FINAL
scalar concaveAngle() const
label nGrow() const
If points get not extruded do nGrow layers of connected faces.
static scalar layerThickness(const thicknessModelType, const label nLayers, const scalar firstLayerThickness, const scalar finalLayerThickness, const scalar totalThickness, const scalar expansionRatio)
Determine overall thickness. Uses two of the four parameters.
label nOuterIter() const
Outer loop to add layer by layer. Can be set to >= max layers.
label nBufferCellsNoExtrude() const
Create buffer region for new layer terminations.
label nLayerIter() const
Number of overall layer addition iterations.
static scalar finalLayerThicknessRatio(const label nLayers, const scalar expansionRatio)
Determine ratio of final layer thickness to.
const dictionary & dict() const
const labelList & numLayers() const
How many layers to add.
static labelPairList findDuplicateFacePairs(const polyMesh &)
Helper routine to find all baffles (two boundary faces.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
label nOldCells() const
Number of old cells.
const labelList & faceMap() const
Old face map.
const pointField & preMotionPoints() const
Pre-motion point positions.
const labelList & cellMap() const
Old cell map.
label nOldFaces() const
Number of old faces.
const labelList & pointMap() const
Old point map.
bool hasMotionPoints() const
Has valid preMotionPoints?
Helper class which maintains intersections of (changing) mesh with (static) surfaces.
static List< labelPair > subsetBaffles(const polyMesh &mesh, const labelList &zoneIDs, const List< labelPair > &baffles)
Subset baffles according to zones.
writeType
Enumeration for what to write. Used as a bit-pattern.
debugType
Enumeration for what to debug. Used as a bit-pattern.
FaceMergeType
Enumeration for what to do with co-planar patch faces on a single.
static autoPtr< indirectPrimitivePatch > makePatch(const polyMesh &, const labelList &)
Create patch from set of patches.
static void updateList(const labelList &newToOld, const T &nullValue, List< T > &elems)
Helper: reorder list according to map.
static writeType writeLevel()
Get/set write level.
OSstream & stream(OSstream *alternative=nullptr)
static void setDisplacement(const labelList &patchIDs, const indirectPrimitivePatch &pp, pointField &patchDisp, pointVectorField &displacement)
Set displacement field from displacement on patch points.
static const complex rootMax
complex (ROOTVGREAT, ROOTVGREAT)
static const complex max
complex (VGREAT,VGREAT)
const word & name() const noexcept
The patch name.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
labelRange range() const noexcept
The face range for all boundary faces.
const labelList & patchID() const
Per boundary face label the patch index.
Mesh consisting of general polyhedral cells.
const faceZoneMesh & faceZones() const noexcept
Return face zone mesh.
void setInstance(const fileName &instance, const IOobject::writeOption wOpt=IOobject::AUTO_WRITE)
Set the instance for mesh files.
virtual const faceList & faces() const
Return raw faces.
bool moving() const noexcept
Is mesh moving.
virtual const labelList & faceOwner() const
Return face owner.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
virtual const pointField & points() const
Return raw points.
A patch is a list of labels that address the faces in the global face list.
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
const labelUList & faceCells() const
Return face-cell addressing.
Direct mesh changes based on v1.3 polyTopoChange syntax.
autoPtr< mapPolyMesh > changeMesh(polyMesh &mesh, const labelUList &patchMap, const bool inflate, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Inplace changes mesh without change of patches.
autoPtr< mapPolyMesh > makeMesh(autoPtr< Type > &newMesh, const IOobject &io, const polyMesh &mesh, const labelUList &patchMap, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Create new mesh with old mesh patches. Additional dictionaries.
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
const labelListList & pointEdges() const
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const vectorField & faceCentres() const
label nInternalFaces() const noexcept
Number of internal faces.
label nPoints() const noexcept
Number of mesh points.
label nCells() const noexcept
Number of mesh cells.
label nFaces() const noexcept
Number of mesh faces.
const cellList & cells() const
label nEdges() const
Number of mesh edges.
int myProcNo() const noexcept
Return processor number.
static word newName(const label myProcNo, const label neighbProcNo)
Return the name of a processorPolyPatch.
All to do with adding layers.
void addLayers(const layerParameters &layerParams, const label nLayerIter, const dictionary &motionDict, const label nRelaxedIter, const label nAllowableErrors, const labelList &patchIDs, const labelList &internalFaceZones, const List< labelPair > &baffles, const labelList &numLayers, const label nIdealTotAddedCells, const globalIndex &globalFaces, indirectPrimitivePatch &pp, const labelListList &edgeGlobalFaces, const labelList &edgePatchID, const labelList &edgeZoneID, const boolList &edgeFlip, const labelList &inflateFaceID, const scalarField &thickness, const scalarIOField &minThickness, const scalarField &expansionRatio, vectorField &patchDisp, labelList &patchNLayers, List< extrudeMode > &extrudeStatus, polyTopoChange &savedMeshMod, labelList &cellNLayers, scalarField &faceRealThickness)
void mergePatchFacesUndo(const layerParameters &layerParams, const dictionary &motionDict, const meshRefinement::FaceMergeType mergeType)
Merge patch faces on same cell.
void doLayers(const dictionary &shrinkDict, const dictionary &motionDict, const layerParameters &layerParams, const meshRefinement::FaceMergeType mergeType, const bool preBalance, decompositionMethod &decomposer, fvMeshDistribute &distributor)
Add layers according to the dictionary settings.
faceZoneType
What to do with faceZone faces.
A class for managing temporary objects.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
const word & name() const noexcept
The zone name.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
const volScalarField & p0
const polyBoundaryMesh & patches
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
surfaceVectorField faceNormals(mesh.Sf()/mesh.magSf())
#define IOWarningInFunction(ios)
Report an IO warning using Foam::Warning.
#define WarningInFunction
Report a warning using Foam::Warning.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
Pair< label > labelPair
A pair of labels.
List< word > wordList
A List of words.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Smanip< ios_base::fmtflags > setf(const ios_base::fmtflags flags)
List< labelPair > labelPairList
List of labelPairs.
Type gSum(const FieldField< Field, Type > &f)
const dimensionSet dimless
Dimensionless.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
List< label > labelList
A List of labels.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
line< point, const point & > linePointRef
A line using referred points.
prefixOSstream Perr
OSstream wrapped stderr (std::cerr) with parallel prefix.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
messageStream Info
Information stream (stdout output on master, null elsewhere)
IOField< scalar > scalarIOField
scalarField with IO.
GeometricField< vector, pointPatchField, pointMesh > pointVectorField
vectorField pointField
pointField is a vectorField.
vector point
Point is a vector.
Omanip< int > setprecision(const int i)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Omanip< int > setw(const int i)
label mergePoints(const PointList &points, labelList &pointToUnique, labelList &uniquePoints, const scalar mergeTol=SMALL, const bool verbose=false)
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Ostream & indent(Ostream &os)
Indent stream.
constexpr scalar degToRad() noexcept
Multiplication factor for degrees to radians conversion.
void inplaceReorder(const labelUList &oldToNew, ListType &input, const bool prune=false)
Inplace reorder the elements of a list.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
List< labelList > labelListList
A List of labelList.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
pointField vertices(const blockVertexList &bvl)
errorManip< error > abort(error &err)
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Field< vector > vectorField
Specialisation of Field<T> for vector.
List< bool > boolList
A List of bools.
static constexpr const zero Zero
Global zero (0)
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
labelListList invertOneToMany(const label len, const labelUList &map)
Invert one-to-many map. Unmapped elements will be size 0.
List< face > faceList
A List of faces.
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
A PrimitivePatch with an IndirectList for the faces, const reference for the point field.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Type gMax(const FieldField< Field, Type > &f)
UIndirectList< label > labelUIndList
UIndirectList of labels.
dimensionedScalar cos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
#define addProfiling(name, descr)
Define profiling trigger with specified name and description string.
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells].reset(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
#define forAll(list, i)
Loop across all elements in list.
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
static const char *const typeName
The type name used in ensight case files.
Unit conversion functions.