65 Foam::label Foam::meshRefinement::createBaffle
70 polyTopoChange& meshMod
73 const face&
f = mesh_.
faces()[faceI];
75 bool zoneFlip =
false;
79 const faceZone& fZone = mesh_.
faceZones()[zoneID];
80 zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
107 <<
"No neighbour patch for internal face " << faceI
112 bool reverseFlip =
false;
115 reverseFlip = !zoneFlip;
118 dupFaceI = meshMod.setAction
178 void Foam::meshRefinement::getIntersections
188 autoPtr<OBJstream> str;
189 if (
debug&OBJINTERSECTIONS)
196 mesh_.time().path()/
timeName()/
"intersections.obj" 200 Pout<<
"getIntersections : Writing surface intersections to file " 205 globalRegion1.setSize(mesh_.nFaces());
207 globalRegion2.setSize(mesh_.nFaces());
235 List<pointIndexHit> hit1;
238 List<pointIndexHit> hit2;
240 surfaces_.findNearestIntersection
258 label faceI = testFaces[i];
260 if (hit1[i].hit() && hit2[i].hit())
264 str().writeLine(start[i], hit1[i].
point());
265 str().writeLine(hit1[i].
point(), hit2[i].
point());
266 str().writeLine(hit2[i].
point(),
end[i]);
270 globalRegion1[faceI] =
271 surfaces_.globalRegion(surface1[i], region1[i]);
272 globalRegion2[faceI] =
273 surfaces_.globalRegion(surface2[i], region2[i]);
275 if (globalRegion1[faceI] == -1 || globalRegion2[faceI] == -1)
285 void Foam::meshRefinement::getBafflePatches
287 const label nErodeCellZones,
292 const bool exitIfLeakPath,
293 const refPtr<coordSetWriter>& leakPathFormatter,
294 refPtr<surfaceWriter>& surfFormatter,
319 bitSet posOrientation;
327 locationsOutsideMesh,
350 ownPatch.setSize(mesh_.nFaces());
352 neiPatch.setSize(mesh_.nFaces());
357 if (unnamedRegion1[faceI] != -1 || unnamedRegion2[faceI] != -1)
359 label ownMasterPatch = -1;
360 if (unnamedRegion1[faceI] != -1)
362 ownMasterPatch = globalToMasterPatch[unnamedRegion1[faceI]];
364 label neiMasterPatch = -1;
365 if (unnamedRegion2[faceI] != -1)
367 neiMasterPatch = globalToMasterPatch[unnamedRegion2[faceI]];
374 label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
377 if (mesh_.isInternalFace(faceI))
379 neiZone = cellToZone[mesh_.faceNeighbour()[faceI]];
383 neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
391 (ownZone >= 0 && neiZone != -2)
392 || (neiZone >= 0 && ownZone != -2)
395 namedSurfaceRegion.size() == 0
396 || namedSurfaceRegion[faceI] == -1
407 ownPatch[faceI] = ownMasterPatch;
408 neiPatch[faceI] = neiMasterPatch;
428 const bool allowBoundary,
433 Map<labelPair> bafflePatch(mesh_.nFaces()/1000);
435 const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();
440 const wordList& faceZoneNames = surfZones[surfI].faceZoneNames();
442 forAll(faceZoneNames, fzi)
445 const word& faceZoneName = faceZoneNames[fzi];
446 label zoneI = fZones.findZoneID(faceZoneName);
447 const faceZone& fZone = fZones[zoneI];
450 label globalRegionI = surfaces_.globalRegion(surfI, fzi);
453 globalToMasterPatch[globalRegionI],
454 globalToSlavePatch[globalRegionI]
457 Info<<
"For zone " << fZone.name() <<
" found patches " 458 << mesh_.boundaryMesh()[zPatches[0]].name() <<
" and " 459 << mesh_.boundaryMesh()[zPatches[1]].name()
464 label faceI = fZone[i];
466 if (allowBoundary || mesh_.isInternalFace(faceI))
469 if (fZone.flipMap()[i])
474 if (!bafflePatch.insert(faceI,
patches))
478 <<
" fc:" << mesh_.faceCentres()[faceI]
479 <<
" in zone " << fZone.name()
480 <<
" is in multiple zones!" 499 ownPatch.size() != mesh_.nFaces()
500 || neiPatch.size() != mesh_.nFaces()
505 <<
" ownPatch:" << ownPatch.size()
506 <<
" neiPatch:" << neiPatch.size()
507 <<
". Should be number of faces:" << mesh_.nFaces()
518 forAll(syncedOwnPatch, faceI)
522 (ownPatch[faceI] == -1 && syncedOwnPatch[faceI] != -1)
523 || (neiPatch[faceI] == -1 && syncedNeiPatch[faceI] != -1)
527 <<
"Non synchronised at face:" << faceI
528 <<
" on patch:" << mesh_.boundaryMesh().whichPatch(faceI)
529 <<
" fc:" << mesh_.faceCentres()[faceI] <<
endl 530 <<
"ownPatch:" << ownPatch[faceI]
531 <<
" syncedOwnPatch:" << syncedOwnPatch[faceI]
532 <<
" neiPatch:" << neiPatch[faceI]
533 <<
" syncedNeiPatch:" << syncedNeiPatch[faceI]
540 polyTopoChange meshMod(mesh_);
544 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
546 if (ownPatch[faceI] != -1)
560 const polyBoundaryMesh&
pbm = mesh_.boundaryMesh();
564 const polyPatch&
pp =
pbm[patchI];
566 label coupledPatchI = -1;
570 && !refCast<const coupledPolyPatch>(
pp).separated()
573 coupledPatchI = patchI;
578 label faceI =
pp.start()+i;
580 if (ownPatch[faceI] != -1)
590 if (coupledPatchI != -1)
592 faceToCoupledPatch_.insert(faceI, coupledPatchI);
601 autoPtr<mapPolyMesh> mapPtr;
609 mapPtr = meshMod.changeMesh(mesh_,
false,
true);
610 mapPolyMesh& map = *mapPtr;
613 mesh_.updateMesh(map);
616 if (map.hasMotionPoints())
618 mesh_.movePoints(map.preMotionPoints());
632 faceSet baffledFacesSet(mesh_,
"baffledFacesSet", 2*nBaffles);
634 const labelList& reverseFaceMap = map.reverseFaceMap();
638 forAll(ownPatch, oldFaceI)
640 label faceI = reverseFaceMap[oldFaceI];
642 if (ownPatch[oldFaceI] != -1 && faceI >= 0)
644 const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
646 baffledFacesSet.insert(ownFaces);
652 label oldFaceI =
faceMap[faceI];
654 if (oldFaceI >= 0 && reverseFaceMap[oldFaceI] != faceI)
656 const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
658 baffledFacesSet.insert(ownFaces);
661 baffledFacesSet.sync(mesh_);
663 updateMesh(map, baffledFacesSet.toc());
681 const faceZone& fZone = faceZones[zoneI];
685 bool hasInfo = getFaceZoneInfo(fZone.
name(), mpI, spI, fzType);
687 if (hasInfo && fzTypes.
found(fzType))
709 for (
const label zoneID :
zoneIDs)
720 if (faceToZone[
p[0]] != -1 && (faceToZone[
p[0]] == faceToZone[
p[1]]))
742 ownPatch.
setSize(mesh_.nFaces());
744 neiPatch.
setSize(mesh_.nFaces());
753 const bitSet isInternalOrCoupled
761 const faceZone& fz = faceZones[zoneI];
762 const word& masterName = faceZoneToMasterPatch_[fz.
name()];
763 label masterPatchI = mesh_.boundaryMesh().findPatchID(masterName);
764 const word& slaveName = faceZoneToSlavePatch_[fz.
name()];
765 label slavePatchI = mesh_.boundaryMesh().findPatchID(slaveName);
767 if (masterPatchI == -1 || slavePatchI == -1)
770 <<
"Problem: masterPatchI:" << masterPatchI
777 if (isInternalOrCoupled[faceI])
781 ownPatch[faceI] = slavePatchI;
782 neiPatch[faceI] = masterPatchI;
786 ownPatch[faceI] = masterPatchI;
787 neiPatch[faceI] = slavePatchI;
812 Info<<
"Converting zoned faces into baffles ..." <<
endl;
821 label nLocalBaffles =
sum(nBaffles);
825 if (nTotalBaffles > 0)
830 <<
setf(ios_base::left)
831 <<
setw(30) <<
"FaceZone" 832 <<
setw(10) <<
"FaceType" 833 <<
setw(10) <<
"nBaffles" 835 <<
setw(30) <<
"--------" 836 <<
setw(10) <<
"--------" 837 <<
setw(10) <<
"--------" 843 const faceZone& fz = faceZones[zoneI];
847 bool hasInfo = getFaceZoneInfo(fz.
name(), mpI, spI, fzType);
855 <<
setw(10) << nBaffles[j]
862 map = createBaffles(ownPatch, neiPatch);
868 baffles.
setSize(nLocalBaffles);
869 originatingFaceZone.
setSize(nLocalBaffles);
873 const labelList& reverseFaceMap = map().reverseFaceMap();
877 label faceI = mesh_.nInternalFaces();
878 faceI < mesh_.nFaces();
882 label oldFaceI =
faceMap[faceI];
883 label masterFaceI = reverseFaceMap[oldFaceI];
884 if (masterFaceI != faceI && ownPatch[oldFaceI] != -1)
886 baffles[baffleI] =
labelPair(masterFaceI, faceI);
887 originatingFaceZone[baffleI] =
faceZoneID[oldFaceI];
892 if (baffleI != baffles.
size())
895 <<
"Had " << baffles.
size() <<
" baffles to create " 896 <<
" but encountered " << baffleI
897 <<
" slave faces originating from patcheable faces." 903 const_cast<Time&
>(mesh_.time())++;
904 Pout<<
"Writing zone-baffled mesh to time " <<
timeName()
910 mesh_.time().path()/
"baffles" 914 Info<<
"Created " << nTotalBaffles <<
" baffles in = " 915 << mesh_.time().cpuTimeIncrement() <<
" s\n" <<
nl <<
endl;
920 originatingFaceZone.
clear();
930 const scalar planarAngle
957 const label baffleValue = 1000000;
972 label faceI =
pp.start();
976 const labelList& fEdges = mesh_.faceEdges(faceI);
980 nBafflesPerEdge[fEdges[fEdgeI]]++;
988 DynamicList<label> fe0;
989 DynamicList<label> fe1;
998 label f0 = couples[i].
first();
999 const labelList& fEdges0 = mesh_.faceEdges(f0, fe0);
1002 nBafflesPerEdge[fEdges0[fEdgeI]] += baffleValue;
1007 label f1 = couples[i].second();
1008 const labelList& fEdges1 = mesh_.faceEdges(f1, fe1);
1011 nBafflesPerEdge[fEdges1[fEdgeI]] += baffleValue;
1030 List<labelPair> filteredCouples(couples.
size());
1043 const labelList& fEdges = mesh_.faceEdges(couple.first());
1047 label edgeI = fEdges[fEdgeI];
1049 if (nBafflesPerEdge[edgeI] == 2*baffleValue+2*1)
1051 filteredCouples[filterI++] = couple;
1057 filteredCouples.setSize(filterI);
1060 label nFiltered =
returnReduce(filteredCouples.size(), sumOp<label>());
1062 Info<<
"freeStandingBaffles : detected " 1064 <<
" free-standing baffles out of " 1077 const pointField& cellCentres = mesh_.cellCentres();
1079 forAll(filteredCouples, i)
1081 const labelPair& couple = filteredCouples[i];
1082 start[i] = cellCentres[mesh_.faceOwner()[couple.first()]];
1083 end[i] = cellCentres[mesh_.faceOwner()[couple.second()]];
1097 List<pointIndexHit> hit1;
1102 List<pointIndexHit> hit2;
1106 surfaces_.findNearestIntersection
1108 identity(surfaces_.surfaces().size()),
1132 forAll(filteredCouples, i)
1134 const labelPair& couple = filteredCouples[i];
1148 && hit1[i].
point().dist(hit2[i].
point()) > mergeDistance_
1159 if ((normal1[i]&normal2[i]) > planarAngleCos)
1163 scalar magN =
mag(
n);
1166 filteredCouples[filterI++] = couple;
1170 else if (hit1[i].hit() || hit2[i].hit())
1176 filteredCouples.setSize(filterI);
1178 Info<<
"freeStandingBaffles : detected " 1180 <<
" planar (within " << planarAngle
1181 <<
" degrees) free-standing baffles out of " 1186 return filteredCouples;
1203 const faceList& faces = mesh_.faces();
1204 const labelList& faceOwner = mesh_.faceOwner();
1209 label face0 = couples[i].
first();
1210 label face1 = couples[i].second();
1215 label own0 = faceOwner[face0];
1216 label own1 = faceOwner[face1];
1218 if (face1 < 0 || own0 < own1)
1221 label zoneID = faceZones.
whichZone(face0);
1222 bool zoneFlip =
false;
1226 const faceZone& fZone = faceZones[zoneID];
1230 label nei = (face1 < 0 ? -1 : own1);
1252 label zoneID = faceZones.
whichZone(face1);
1253 bool zoneFlip =
false;
1257 const faceZone& fZone = faceZones[zoneID];
1282 const label faceI = iter.key();
1283 const label patchI = iter.val();
1285 if (!mesh_.isInternalFace(faceI))
1288 <<
"problem: face:" << faceI
1289 <<
" at:" << mesh_.faceCentres()[faceI]
1290 <<
"(wanted patch:" << patchI
1294 label zoneID = faceZones.
whichZone(faceI);
1295 bool zoneFlip =
false;
1299 const faceZone& fZone = faceZones[zoneID];
1323 mesh_.moving(
false);
1326 mapPtr = meshMod.
changeMesh(mesh_,
false,
true);
1330 mesh_.updateMesh(map);
1357 newExposedFaces[newI++] = newFace0;
1360 const label newFace1 = map.
reverseFaceMap()[couples[i].second()];
1363 newExposedFaces[newI++] = newFace1;
1366 newExposedFaces.
setSize(newI);
1367 updateMesh(map, newExposedFaces);
1376 const bool doInternalZones,
1377 const bool doBaffleZones
1383 if (doInternalZones)
1407 mapPtr = mergeBaffles(zoneBaffles,
Map<label>(0));
1414 void Foam::meshRefinement::findCellZoneGeometric
1424 const pointField& cellCentres = mesh_.cellCentres();
1425 const labelList& faceOwner = mesh_.faceOwner();
1426 const labelList& faceNeighbour = mesh_.faceNeighbour();
1430 surfaces_.findInside
1432 closedNamedSurfaces,
1437 forAll(insideSurfaces, cellI)
1439 label surfI = insideSurfaces[cellI];
1443 if (cellToZone[cellI] == -2)
1445 cellToZone[cellI] = surfaceToCellZone[surfI];
1447 else if (cellToZone[cellI] == -1)
1452 cellToZone[cellI] = surfaceToCellZone[surfI];
1465 label nCandidates = 0;
1466 forAll(namedSurfaceRegion, faceI)
1468 if (namedSurfaceRegion[faceI] != -1)
1470 if (mesh_.isInternalFace(faceI))
1484 forAll(namedSurfaceRegion, faceI)
1486 if (namedSurfaceRegion[faceI] != -1)
1488 label own = faceOwner[faceI];
1489 const point& ownCc = cellCentres[own];
1491 if (mesh_.isInternalFace(faceI))
1493 label nei = faceNeighbour[faceI];
1494 const point& neiCc = cellCentres[nei];
1496 const vector d = 1
e-4*(neiCc - ownCc);
1497 candidatePoints[nCandidates++] = ownCc-d;
1498 candidatePoints[nCandidates++] = neiCc+d;
1503 const point& neiFc = neiCc[faceI-mesh_.nInternalFaces()];
1506 const vector d = 1
e-4*(neiFc - ownCc);
1507 candidatePoints[nCandidates++] = ownCc-d;
1515 surfaces_.findInside
1517 closedNamedSurfaces,
1526 forAll(namedSurfaceRegion, faceI)
1528 if (namedSurfaceRegion[faceI] != -1)
1530 label own = faceOwner[faceI];
1532 if (mesh_.isInternalFace(faceI))
1534 label ownSurfI = insideSurfaces[nCandidates++];
1537 cellToZone[own] = surfaceToCellZone[ownSurfI];
1540 label neiSurfI = insideSurfaces[nCandidates++];
1543 label nei = faceNeighbour[faceI];
1545 cellToZone[nei] = surfaceToCellZone[neiSurfI];
1550 label ownSurfI = insideSurfaces[nCandidates++];
1553 cellToZone[own] = surfaceToCellZone[ownSurfI];
1564 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1566 label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
1567 label neiZone = cellToZone[mesh_.faceNeighbour()[faceI]];
1569 if (namedSurfaceRegion[faceI] == -1 && (ownZone != neiZone))
1579 else if (neiZone == -1)
1585 minZone =
min(ownZone, neiZone);
1590 label geomSurfI = surfaceToCellZone.
find(minZone);
1592 if (geomSurfI != -1)
1594 namedSurfaceRegion[faceI] =
1595 surfaces_.globalRegion(geomSurfI, 0);
1603 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
1613 label faceI =
pp.start()+i;
1614 label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
1615 label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
1617 if (namedSurfaceRegion[faceI] == -1 && (ownZone != neiZone))
1625 else if (neiZone == -1)
1631 minZone =
min(ownZone, neiZone);
1636 label geomSurfI = surfaceToCellZone.
find(minZone);
1638 if (geomSurfI != -1)
1640 namedSurfaceRegion[faceI] =
1641 surfaces_.globalRegion(geomSurfI, 0);
1653 void Foam::meshRefinement::findCellZoneInsideWalk
1663 boolList blockedFace(mesh_.nFaces());
1666 forAll(blockedFace, faceI)
1668 if (faceToZone[faceI] == -1)
1670 blockedFace[faceI] =
false;
1674 blockedFace[faceI] =
true;
1680 regionSplit cellRegion(mesh_, blockedFace);
1681 blockedFace.clear();
1684 (void)mesh_.tetBasePtIs();
1687 forAll(locationsInMesh, i)
1690 const point& insidePoint = locationsInMesh[i];
1691 label zoneID = zonesInMesh[i];
1694 label keepRegionI = findRegion
1698 vector::uniform(mergeDistance_),
1702 Info<<
"For cellZone " 1706 : mesh_.cellZones()[zoneID].name()
1708 <<
" found point " << insidePoint
1709 <<
" in global region " << keepRegionI
1710 <<
" out of " << cellRegion.nRegions() <<
" regions." <<
endl;
1712 if (keepRegionI == -1)
1715 <<
"Point " << insidePoint
1716 <<
" is not inside the mesh." <<
nl 1717 <<
"Bounding box of the mesh:" << mesh_.bounds()
1726 label nWarnings = 0;
1728 forAll(cellRegion, cellI)
1730 if (cellRegion[cellI] == keepRegionI)
1732 if (cellToZone[cellI] == -2)
1735 cellToZone[cellI] = zoneID;
1737 else if (cellToZone[cellI] != zoneID)
1739 if (cellToZone[cellI] != -1 && nWarnings < 10)
1743 <<
" at " << mesh_.cellCentres()[cellI]
1744 <<
" is inside cellZone " 1748 : mesh_.cellZones()[zoneID].name()
1750 <<
" from locationInMesh " << insidePoint
1751 <<
" but already marked as being in zone " 1752 << mesh_.cellZones()[cellToZone[cellI]].name()
1754 <<
"This can happen if your surfaces are not" 1755 <<
" (sufficiently) closed." 1761 cellToZone[cellI] = zoneID;
1769 void Foam::meshRefinement::findCellZoneInsideWalk
1781 forAll(zoneNamesInMesh, i)
1783 zoneIDs[i] = czs.findZoneID(zoneNamesInMesh[i]);
1785 findCellZoneInsideWalk
1795 bool Foam::meshRefinement::calcRegionToZone
1797 const label backgroundZoneID,
1798 const label surfZoneI,
1799 const label ownRegion,
1800 const label neiRegion,
1805 bool changed =
false;
1808 if (ownRegion != neiRegion)
1815 if (regionToCellZone[ownRegion] == -2)
1817 if (surfZoneI == -1)
1822 if (regionToCellZone[neiRegion] != -2)
1824 regionToCellZone[ownRegion] = regionToCellZone[neiRegion];
1828 else if (regionToCellZone[neiRegion] == surfZoneI)
1833 if (backgroundZoneID != -2)
1835 regionToCellZone[ownRegion] = backgroundZoneID;
1839 else if (regionToCellZone[neiRegion] != -2)
1843 regionToCellZone[ownRegion] = surfZoneI;
1847 else if (regionToCellZone[neiRegion] == -2)
1849 if (surfZoneI == -1)
1854 regionToCellZone[neiRegion] = regionToCellZone[ownRegion];
1857 else if (regionToCellZone[ownRegion] == surfZoneI)
1861 if (backgroundZoneID != -2)
1863 regionToCellZone[neiRegion] = backgroundZoneID;
1867 else if (regionToCellZone[ownRegion] != -2)
1871 regionToCellZone[neiRegion] = surfZoneI;
1880 void Foam::meshRefinement::findCellZoneTopo
1882 const label backgroundZoneID,
1911 boolList blockedFace(mesh_.nFaces());
1913 forAll(unnamedSurfaceRegion, faceI)
1917 unnamedSurfaceRegion[faceI] == -1
1918 && namedSurfaceRegion[faceI] == -1
1921 blockedFace[faceI] =
false;
1925 blockedFace[faceI] =
true;
1931 regionSplit cellRegion(mesh_, blockedFace);
1932 blockedFace.clear();
1938 labelList regionToCellZone(cellRegion.nRegions(), -2);
1943 forAll(cellToZone, cellI)
1945 label regionI = cellRegion[cellI];
1946 if (cellToZone[cellI] != -2 && regionToCellZone[regionI] == -2)
1948 regionToCellZone[regionI] = cellToZone[cellI];
1961 forAll(locationsInMesh, i)
1963 const point& keepPoint = locationsInMesh[i];
1964 label keepRegionI = findRegion
1968 vector::uniform(mergeDistance_),
1972 Info<<
"Found point " << keepPoint
1973 <<
" in global region " << keepRegionI
1974 <<
" out of " << cellRegion.nRegions() <<
" regions." <<
endl;
1976 if (keepRegionI == -1)
1979 <<
"Point " << keepPoint
1980 <<
" is not inside the mesh." <<
nl 1981 <<
"Bounding box of the mesh:" << mesh_.bounds()
1989 if (regionToCellZone[keepRegionI] == -2)
1991 regionToCellZone[keepRegionI] = -1;
2010 bool changed =
false;
2014 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
2016 label regionI = namedSurfaceRegion[faceI];
2019 if (unnamedSurfaceRegion[faceI] == -1 && regionI != -1)
2024 label surfI = surfaces_.whichSurface(regionI).first();
2026 bool changedCell = calcRegionToZone
2029 surfaceToCellZone[surfI],
2030 cellRegion[mesh_.faceOwner()[faceI]],
2031 cellRegion[mesh_.faceNeighbour()[faceI]],
2035 changed = changed | changedCell;
2041 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
2057 label faceI =
pp.start()+i;
2059 label regionI = namedSurfaceRegion[faceI];
2062 if (unnamedSurfaceRegion[faceI] == -1 && regionI != -1)
2064 label surfI = surfaces_.whichSurface(regionI).first();
2066 bool changedCell = calcRegionToZone
2069 surfaceToCellZone[surfI],
2070 cellRegion[mesh_.faceOwner()[faceI]],
2071 neiCellRegion[faceI-mesh_.nInternalFaces()],
2075 changed = changed | changedCell;
2090 Pout<<
"meshRefinement::findCellZoneTopo :" 2091 <<
" nRegions:" << regionToCellZone.size()
2092 <<
" of which visited (-1 = background, >= 0 : cellZone) :" 2095 forAll(regionToCellZone, regionI)
2097 if (regionToCellZone[regionI] != -2)
2099 Pout<<
"Region " << regionI
2100 <<
" becomes cellZone:" << regionToCellZone[regionI]
2107 forAll(cellToZone, cellI)
2109 label regionI = cellRegion[cellI];
2110 if (cellToZone[cellI] == -2 && regionToCellZone[regionI] != -2)
2112 cellToZone[cellI] = regionToCellZone[regionI];
2118 void Foam::meshRefinement::erodeCellZone
2120 const label nErodeCellZones,
2121 const label backgroundZoneID,
2138 for (label iter = 0; iter < nErodeCellZones; iter++)
2145 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
2149 unnamedSurfaceRegion[facei] == -1
2150 && namedSurfaceRegion[facei] == -1
2153 label own = mesh_.faceOwner()[facei];
2154 label nei = mesh_.faceNeighbour()[facei];
2155 if (cellToZone[own] == -2 && cellToZone[nei] >= -1)
2157 erodedCellToZone[nei] = backgroundZoneID;
2160 else if (cellToZone[nei] == -2 && cellToZone[own] >= -1)
2162 erodedCellToZone[own] = backgroundZoneID;
2170 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
2186 label facei =
pp.start()+i;
2189 unnamedSurfaceRegion[facei] == -1
2190 && namedSurfaceRegion[facei] == -1
2193 label own = mesh_.faceOwner()[facei];
2194 label bFacei = facei-mesh_.nInternalFaces();
2195 if (neiCellZone[bFacei] == -2 && cellToZone[own] >= -1)
2197 erodedCellToZone[own] = backgroundZoneID;
2205 cellToZone.transfer(erodedCellToZone);
2207 reduce(nChanged, sumOp<label>());
2210 Pout<<
"erodeCellZone : eroded " << nChanged
2211 <<
" cells (moved from cellZone to background zone " 2212 << backgroundZoneID <<
endl;
2223 void Foam::meshRefinement::growCellZone
2225 const label nGrowCellZones,
2226 const label backgroundZoneID,
2233 if (nGrowCellZones != 1)
2236 <<
"Growing currently only supported for single layer." 2237 <<
" Exiting zone growing" <<
endl;
2254 const FixedList<label, 3> wallTag
2264 List<FixedList<label, 3>> surfaces(1);
2285 forAll(cellToZone, celli)
2287 if (cellToZone[celli] >= 0)
2289 const FixedList<label, 3> zoneTag
2296 origins[0] = mesh_.cellCentres()[celli];
2298 surfaces[0] = zoneTag;
2299 allCellInfo[celli] = wallPoints(origins, distSqrs, surfaces);
2304 DynamicList<wallPoints> faceDist(mesh_.nFaces()/128);
2305 DynamicList<label> changedFaces(mesh_.nFaces()/128);
2307 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
2309 const label own = mesh_.faceOwner()[facei];
2310 const label nei = mesh_.faceNeighbour()[facei];
2311 if (cellToZone[own] >= 0 && cellToZone[nei] < 0)
2315 origins[0] = mesh_.faceCentres()[facei];
2317 surfaces[0] = FixedList<label, 3>
2323 faceDist.append(wallPoints(origins, distSqrs, surfaces));
2324 changedFaces.append(facei);
2326 else if (cellToZone[own] < 0 && cellToZone[nei] >= 0)
2330 origins[0] = mesh_.faceCentres()[facei];
2332 surfaces[0] = FixedList<label, 3>
2338 faceDist.append(wallPoints(origins, distSqrs, surfaces));
2339 changedFaces.append(facei);
2343 unnamedSurfaceRegion1[facei] != -1
2344 || unnamedSurfaceRegion2[facei] != -1
2349 origins[0] = mesh_.faceCentres()[facei];
2351 surfaces[0] = wallTag;
2352 faceDist.append(wallPoints(origins, distSqrs, surfaces));
2353 changedFaces.append(facei);
2361 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
2374 label facei =
pp.start()+i;
2375 label own = mesh_.faceOwner()[facei];
2376 label bFacei = facei-mesh_.nInternalFaces();
2377 if (cellToZone[own] >= 0 && neiCellZone[bFacei] < 0)
2381 surfaces[0] = FixedList<label, 3>
2387 faceDist.append(wallPoints(origins, distSqrs, surfaces));
2388 changedFaces.append(facei);
2390 else if (cellToZone[own] < 0 && neiCellZone[bFacei] >= 0)
2396 unnamedSurfaceRegion1[facei] != -1
2397 || unnamedSurfaceRegion2[facei] != -1
2401 origins[0] = mesh_.faceCentres()[facei];
2403 surfaces[0] = wallTag;
2404 faceDist.append(wallPoints(origins, distSqrs, surfaces));
2405 changedFaces.append(facei);
2414 label facei =
pp.start()+i;
2415 label own = mesh_.faceOwner()[facei];
2416 if (cellToZone[own] < 0)
2420 surfaces[0] = wallTag;
2421 faceDist.append(wallPoints(origins, distSqrs, surfaces));
2422 changedFaces.append(facei);
2433 List<scalarList> regionToBlockSize(surfaces_.surfaces().size());
2435 forAll(surfaces_.surfaces(), surfi)
2437 const label geomi = surfaces_.surfaces()[surfi];
2438 const auto&
s = surfaces_.geometry()[geomi];
2439 const label nRegions =
s.regions().size();
2440 regionToBlockSize[surfi].setSize(nRegions,
Foam::sqr(GREAT));
2446 FaceCellWave<wallPoints, wallPoints::trackData>
wallDistCalc 2463 bitSet isChangedCell(mesh_.nCells());
2469 const List<FixedList<label, 3>>& surfaces =
2472 if (surfaces.size())
2476 isChangedCell.set(celli);
2479 if (surfaces.size() > 1)
2485 for (label i = 0; i < surfaces.size(); i++)
2487 const label zonei = surfaces[i][0];
2488 const scalar d2 =
allCellInfo[celli].distSqr()[i];
2489 if (zonei != cellToZone[celli] && d2 < minDistSqr)
2498 if (minZone != cellToZone[celli] && minZone != wallTag[0])
2500 cellToZone[celli] = minZone;
2501 isChangedCell.set(celli);
2509 if (backgroundZoneID != -2)
2511 forAll(cellToZone, celli)
2513 if (cellToZone[celli] == -2)
2515 cellToZone[celli] = backgroundZoneID;
2528 for (
const label celli : isChangedCell)
2530 const cell& cFaces = mesh_.cells()[celli];
2531 for (
const label facei : cFaces)
2533 const label own = mesh_.faceOwner()[facei];
2534 const label ownZone = cellToZone[own];
2537 if (mesh_.isInternalFace(facei))
2539 const label neiZone = cellToZone[mesh_.faceNeighbour()[facei]];
2540 nbrZone = (own != celli ? ownZone : neiZone);
2544 nbrZone = neiCellZone[facei-mesh_.nInternalFaces()];
2547 if (nbrZone == cellToZone[celli])
2551 unnamedSurfaceRegion1[facei] != -1
2552 || unnamedSurfaceRegion2[facei] != -1
2555 unnamedSurfaceRegion1[facei] = -1;
2556 unnamedSurfaceRegion2[facei] = -1;
2561 namedSurfaceRegion.size()
2562 && namedSurfaceRegion[facei] != -1
2565 namedSurfaceRegion[facei] = -1;
2572 reduce(nUnnamed, sumOp<label>());
2573 reduce(nNamed, sumOp<label>());
2579 unnamedSurfaceRegion1,
2585 unnamedSurfaceRegion2,
2588 if (namedSurfaceRegion.size())
2600 Pout<<
"growCellZone : grown cellZones by " 2602 <<
" cells (moved from background to nearest cellZone)" 2604 Pout<<
"growCellZone : unmarked " << nUnnamed
2605 <<
" unzoned intersections; " << nNamed <<
" zoned intersections; " 2611 void Foam::meshRefinement::makeConsistentFaceIndex
2626 const labelList& faceOwner = mesh_.faceOwner();
2627 const labelList& faceNeighbour = mesh_.faceNeighbour();
2629 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
2631 label ownZone = cellToZone[faceOwner[faceI]];
2632 label neiZone = cellToZone[faceNeighbour[faceI]];
2633 label globalI = namedSurfaceRegion[faceI];
2639 && surfaceMap[surfaces_.whichSurface(globalI).first()] == -1
2642 namedSurfaceRegion[faceI] = -1;
2646 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
2661 label faceI =
pp.start()+i;
2663 label ownZone = cellToZone[faceOwner[faceI]];
2664 label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
2665 label globalI = namedSurfaceRegion[faceI];
2671 && surfaceMap[surfaces_.whichSurface(globalI).first()] == -1
2674 namedSurfaceRegion[faceI] = -1;
2683 label faceI =
pp.start()+i;
2684 label globalI = namedSurfaceRegion[faceI];
2689 && surfaceMap[surfaces_.whichSurface(globalI).first()] == -1
2692 namedSurfaceRegion[faceI] = -1;
2700 void Foam::meshRefinement::getIntersections
2707 bitSet& posOrientation
2710 namedSurfaceRegion.setSize(mesh_.nFaces());
2711 namedSurfaceRegion = -1;
2713 posOrientation.setSize(mesh_.nFaces());
2714 posOrientation =
false;
2745 List<pointIndexHit> hit1;
2749 List<pointIndexHit> hit2;
2752 surfaces_.findNearestIntersection
2771 label faceI = testFaces[i];
2772 const vector&
area = mesh_.faceAreas()[faceI];
2774 if (surface1[i] != -1)
2781 magSqr(hit2[i].hitPoint())
2782 <
magSqr(hit1[i].hitPoint())
2786 namedSurfaceRegion[faceI] = surfaces_.globalRegion
2791 posOrientation.
set(faceI, ((
area&normal2[i]) > 0));
2792 nSurfFaces[surface2[i]]++;
2796 namedSurfaceRegion[faceI] = surfaces_.globalRegion
2801 posOrientation.set(faceI, ((
area&normal1[i]) > 0));
2802 nSurfFaces[surface1[i]]++;
2805 else if (surface2[i] != -1)
2807 namedSurfaceRegion[faceI] = surfaces_.globalRegion
2812 posOrientation.set(faceI, ((
area&normal2[i]) > 0));
2813 nSurfFaces[surface2[i]]++;
2832 forAll(nSurfFaces, surfI)
2835 << surfaces_.names()[surfI]
2836 <<
" nZoneFaces:" << nSurfFaces[surfI] <<
nl;
2843 void Foam::meshRefinement::zonify
2845 const bool allowFreeStandingZoneFaces,
2846 const label nErodeCellZones,
2847 const label backgroundZoneID,
2851 const bool exitIfLeakPath,
2852 const refPtr<coordSetWriter>& leakPathFormatter,
2853 refPtr<surfaceWriter>& surfFormatter,
2859 bitSet& posOrientation
2872 const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();
2875 const List<pointField> allLocations
2881 locationsOutsideMesh
2886 labelList neiLevel(mesh_.nBoundaryFaces());
2888 calcNeighbourData(neiLevel, neiCc);
2896 if (namedSurfaces.size())
2908 cellToZone.setSize(mesh_.nCells());
2911 namedSurfaceRegion.clear();
2912 posOrientation.clear();
2932 autoPtr<mapDistribute> unnamedMapPtr;
2933 if (locationsOutsideMesh.size())
2938 [](
const label
x){
return x != -1;}
2941 const globalIndex globalUnnamedFaces(unnamedFaces.size());
2951 unnamedClosureFaces,
2957 Pout<<
"meshRefinement::zonify : found wall closure faces:" 2958 << unnamedClosureFaces.size()
2959 <<
" map:" << bool(unnamedMapPtr) <<
endl;
2967 if (leakPathFormatter)
2969 boolList blockedFace(mesh_.nFaces(),
false);
2970 UIndirectList<bool>(blockedFace, unnamedFaces) =
true;
2971 const fileName fName
2977 locationsOutsideMesh,
2979 leakPathFormatter.constCast()
2982 Info<<
"Dumped leak path to " << fName <<
endl;
2992 err <<
"Locations in mesh " << locationsInMesh
2993 <<
" connect to one of the locations outside mesh " 2994 << locationsOutsideMesh <<
endl;
3004 UIndirectList<label>(unnamedRegion1, unnamedFaces)
3006 unnamedMapPtr->distribute(packedRegion1);
3009 UIndirectList<label>(unnamedRegion2, unnamedFaces)
3011 unnamedMapPtr->distribute(packedRegion2);
3012 forAll(unnamedClosureFaces, i)
3014 const label sloti = unnamedToClosure[i];
3018 const label facei = unnamedClosureFaces[i];
3019 const label region1 = unnamedRegion1[facei];
3020 const label slotRegion1 = packedRegion1[sloti];
3021 const label region2 = unnamedRegion2[facei];
3022 const label slotRegion2 = packedRegion2[sloti];
3024 if (slotRegion1 != region1 || slotRegion2 != region2)
3026 unnamedRegion1[facei] = slotRegion1;
3027 unnamedRegion2[facei] = slotRegion2;
3039 autoPtr<mapDistribute> namedMapPtr;
3040 if (namedSurfaces.size())
3169 bitSet isClosureFace(mesh_.nFaces());
3170 isClosureFace.set(unnamedClosureFaces);
3171 isClosureFace.set(namedClosureFaces);
3172 const labelList closureFaces(isClosureFace.sortedToc());
3176 UIndirectList<face>(mesh_.faces(), closureFaces),
3184 bitSet isFrozenPoint(mesh_.nPoints());
3185 forAll(nEdgeFaces, edgei)
3187 if (nEdgeFaces[edgei] != 1)
3199 const_cast<meshRefinement&
>(*this).addPointZone(
"frozenPoints");
3200 const bitSet oldSet(mesh_.nPoints(), pointZones[zonei]);
3201 isFrozenPoint.
set(oldSet);
3207 orEqOp<unsigned int>(),
3212 pointZones.clearAddressing();
3213 pointZones[zonei] = isFrozenPoint.sortedToc();
3217 mkDir(mesh_.time().timePath());
3218 const pointZone& pz = pointZones[zonei];
3219 OBJstream str(mesh_.time().timePath()/pz.name()+
".obj");
3220 Pout<<
"Writing " << pz.size() <<
" frozen points to " 3222 for (
const label pointi : pz)
3224 str.
write(mesh_.points()[pointi]);
3232 mesh_.time().globalPath()
3234 / mesh_.pointsInstance()
3235 /
"unnamedClosureFaces" 3240 <<
returnReduce(unnamedClosureFaces.size(), sumOp<label>())
3245 IndirectList<face>(mesh_.faces(), unnamedClosureFaces),
3251 setPatch.localPoints(),
3252 setPatch.localFaces(),
3255 surfFormatter->write();
3256 surfFormatter->clear();
3262 mesh_.time().globalPath()
3264 / mesh_.pointsInstance()
3265 /
"namedClosureFaces" 3270 <<
returnReduce(namedClosureFaces.size(), sumOp<label>())
3275 IndirectList<face>(mesh_.faces(), namedClosureFaces),
3281 setPatch.localPoints(),
3282 setPatch.localFaces(),
3285 surfFormatter->write();
3286 surfFormatter->clear();
3295 if (locationsInMesh.size())
3297 Info<<
"Setting cellZones according to locationsInMesh:" <<
endl;
3299 labelList locationsZoneIDs(zonesInMesh.size(), -1);
3300 forAll(locationsInMesh, i)
3302 const word&
name = zonesInMesh[i];
3304 Info<<
"Location : " << locationsInMesh[i] <<
nl 3309 label zoneID = mesh_.cellZones().findZoneID(
name);
3314 locationsZoneIDs[i] = zoneID;
3321 findCellZoneInsideWalk
3338 if (locationSurfaces.size())
3340 Info<<
"Found " << locationSurfaces.size()
3341 <<
" named surfaces with a provided inside point." 3342 <<
" Assigning cells inside these surfaces" 3343 <<
" to the corresponding cellZone." 3349 DynamicField<point>
insidePoints(locationSurfaces.size());
3350 DynamicList<label> insidePointCellZoneIDs(locationSurfaces.size());
3351 forAll(locationSurfaces, i)
3353 const label surfI = locationSurfaces[i];
3354 const auto& surfInsidePoints = surfZones[surfI].zoneInsidePoints();
3356 const word&
name = surfZones[surfI].cellZoneName();
3360 zoneID = mesh_.cellZones().findZoneID(
name);
3364 <<
"Specified non-existing cellZone " <<
name 3365 <<
" for surface " << surfaces_.names()[surfI]
3370 for (
const auto& pt : surfInsidePoints)
3373 insidePointCellZoneIDs.append(zoneID);
3379 labelList allRegion1(mesh_.nFaces(), -1);
3380 forAll(namedSurfaceRegion, faceI)
3382 allRegion1[faceI] =
max 3384 unnamedRegion1[faceI],
3385 namedSurfaceRegion[faceI]
3389 findCellZoneInsideWalk
3392 insidePointCellZoneIDs,
3408 surfaces_.geometry(),
3409 surfaces_.surfaces()
3413 if (closedNamedSurfaces.
size())
3415 Info<<
"Found " << closedNamedSurfaces.
size()
3416 <<
" closed, named surfaces. Assigning cells in/outside" 3417 <<
" these surfaces to the corresponding cellZone." 3420 findCellZoneGeometric
3423 closedNamedSurfaces,
3434 if (namedSurfaces.size())
3436 if (nErodeCellZones == 0)
3438 Info<<
"Walking from known cellZones; crossing a faceZone " 3439 <<
"face changes cellZone" <<
nl <<
endl;
3452 else if (nErodeCellZones < 0)
3454 Info<<
"Growing cellZones by " << -nErodeCellZones
3455 <<
" layers" <<
nl <<
endl;
3469 Info<<
"Eroding cellZone cells to make these consistent with" 3470 <<
" faceZone faces" <<
nl <<
endl;
3486 if (!allowFreeStandingZoneFaces)
3488 Info<<
"Only keeping zone faces inbetween different cellZones." 3502 labelList surfaceMap(surfZones.size(), -1);
3503 forAll(standaloneNamedSurfaces, i)
3505 surfaceMap[standaloneNamedSurfaces[i]] = i;
3508 makeConsistentFaceIndex
3516 else if (nErodeCellZones < 0 &&
gMax(cellToZone) >= 0)
3521 Info<<
"Growing cellZones by " << -nErodeCellZones
3522 <<
" layers" <<
nl <<
endl;
3535 if (!allowFreeStandingZoneFaces && namedSurfaceRegion.size())
3537 Info<<
"Only keeping zone faces inbetween different cellZones." 3551 labelList surfaceMap(surfZones.size(), -1);
3552 forAll(standaloneNamedSurfaces, i)
3554 surfaceMap[standaloneNamedSurfaces[i]] = i;
3557 makeConsistentFaceIndex
3570 label nZones =
gMax(cellToZone)+1;
3572 label nUnvisited = 0;
3573 label nBackgroundCells = 0;
3575 forAll(cellToZone, cellI)
3577 label zoneI = cellToZone[cellI];
3580 nZoneCells[zoneI]++;
3582 else if (zoneI == -1)
3586 else if (zoneI == -2)
3596 reduce(nUnvisited, sumOp<label>());
3597 reduce(nBackgroundCells, sumOp<label>());
3598 forAll(nZoneCells, zoneI)
3600 reduce(nZoneCells[zoneI], sumOp<label>());
3602 Info<<
"nUnvisited :" << nUnvisited <<
endl;
3603 Info<<
"nBackgroundCells:" << nBackgroundCells <<
endl;
3604 Info<<
"nZoneCells :" << nZoneCells <<
endl;
3608 const_cast<Time&
>(mesh_.time())++;
3609 Pout<<
"Writing cell zone allocation on mesh to time " 3614 writeType(writeLevel() | WRITEMESH),
3615 mesh_.time().path()/
"cell2Zone" 3633 forAll(cellToZone, cellI)
3635 volCellToZone[cellI] = cellToZone[cellI];
3637 volCellToZone.write();
3739 void Foam::meshRefinement::handleSnapProblems
3741 const snapParameters& snapParams,
3742 const bool useTopologicalSnapDetection,
3743 const bool removeEdgeConnectedCells,
3745 const dictionary& motionDict,
3752 <<
"Introducing baffles to block off problem cells" <<
nl 3753 <<
"----------------------------------------------" <<
nl 3758 if (useTopologicalSnapDetection)
3760 markFacesOnProblemCells
3763 removeEdgeConnectedCells,
3765 globalToMasterPatch,
3773 markFacesOnProblemCellsGeometric
3777 globalToMasterPatch,
3784 Info<<
"Analyzed problem cells in = " 3789 faceSet problemFaces(mesh_,
"problemFaces", mesh_.nFaces()/100);
3793 if (facePatch[faceI] != -1)
3795 problemFaces.insert(faceI);
3798 problemFaces.instance() =
timeName();
3799 Pout<<
"Dumping " << problemFaces.size()
3800 <<
" problem faces to " << problemFaces.objectPath() <<
endl;
3801 problemFaces.
write();
3804 Info<<
"Introducing baffles to delete problem cells." <<
nl <<
endl;
3815 List<DynamicList<label>> zoneToFaces(fzs.size());
3816 List<DynamicList<bool>> zoneToFlip(fzs.size());
3821 zoneToFaces[zonei].append(fzs[zonei]);
3822 zoneToFlip[zonei].append(fzs[zonei].flipMap());
3828 if (facePatch[facei] != -1)
3830 label zonei = faceZone[facei];
3833 zoneToFaces[zonei].append(facei);
3834 zoneToFlip[zonei].append(
false);
3839 forAll(zoneToFaces, zonei)
3853 createBaffles(facePatch, facePatch);
3860 Info<<
"Created baffles in = " 3863 printMeshInfo(
debug,
"After introducing baffles",
true);
3867 const_cast<Time&
>(mesh_.time())++;
3868 Pout<<
"Writing extra baffled mesh to time " 3873 writeType(writeLevel() | WRITEMESH),
3876 Pout<<
"Dumped debug data in = " 3889 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
3890 const labelList& faceOwner = mesh_.faceOwner();
3891 const labelList& faceNeighbour = mesh_.faceNeighbour();
3908 DynamicList<label>
faceLabels(mesh_.nFaces()/100);
3910 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
3912 if (faceToZone[faceI] != -1)
3915 label ownZone = cellToZone[faceOwner[faceI]];
3916 label neiZone = cellToZone[faceNeighbour[faceI]];
3917 if (ownZone == neiZone)
3929 label faceI =
pp.start()+i;
3930 if (faceToZone[faceI] != -1)
3933 label ownZone = cellToZone[faceOwner[faceI]];
3934 label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
3935 if (ownZone == neiZone)
3946 void Foam::meshRefinement::calcPatchNumMasterFaces
3948 const bitSet& isMasterFace,
3954 nMasterFacesPerEdge.setSize(
patch.nEdges());
3955 nMasterFacesPerEdge = 0;
3959 const label meshFaceI =
patch.addressing()[faceI];
3961 if (isMasterFace.test(meshFaceI))
3966 nMasterFacesPerEdge[fEdges[fEdgeI]]++;
3974 patch.meshEdges(mesh_.edges(), mesh_.pointEdges()),
3975 nMasterFacesPerEdge,
3982 Foam::label Foam::meshRefinement::markPatchZones
3989 List<edgeTopoDistanceData<label>> allEdgeInfo(
patch.nEdges());
3997 forAll(nMasterFacesPerEdge, edgeI)
3999 if (nMasterFacesPerEdge[edgeI] > 2)
4001 allEdgeInfo[edgeI] = edgeTopoDistanceData<label>(0, -2);
4013 DynamicList<label> changedEdges;
4014 DynamicList<edgeTopoDistanceData<label>> changedInfo;
4016 const scalar tol = PatchEdgeFaceWave
4019 edgeTopoDistanceData<label>
4020 >::propagationTol();
4024 const globalIndex globalFaces(
patch.size());
4028 label currentZoneI = 0;
4038 globalSeed = globalFaces.toGlobal(faceI);
4043 reduce(globalSeed, minOp<label>());
4050 if (globalFaces.isLocal(globalSeed))
4052 const label seedFaceI = globalFaces.toLocal(globalSeed);
4054 edgeTopoDistanceData<label>& faceInfo =
allFaceInfo[seedFaceI];
4057 faceInfo = edgeTopoDistanceData<label>(0, currentZoneI);
4063 label edgeI = fEdges[fEdgeI];
4065 edgeTopoDistanceData<label>& edgeInfo = allEdgeInfo[edgeI];
4069 edgeInfo.updateEdge<
int>
4081 changedEdges.append(edgeI);
4082 changedInfo.append(edgeInfo);
4098 edgeTopoDistanceData<label>
4114 faceToZone.setSize(
patch.size());
4120 <<
"Problem: unvisited face " << faceI
4121 <<
" at " <<
patch.faceCentres()[faceI]
4127 return currentZoneI;
4131 void Foam::meshRefinement::consistentOrientation
4133 const bitSet& isMasterFace,
4137 const Map<label>& zoneToOrientation,
4141 const polyBoundaryMesh& bm = mesh_.boundaryMesh();
4144 List<patchFaceOrientation> allEdgeInfo(
patch.nEdges());
4155 const label meshFaceI =
patch.addressing()[faceI];
4156 const label patchI = bm.whichPatch(meshFaceI);
4162 && !isMasterFace.test(meshFaceI)
4175 label nProtected = 0;
4177 forAll(nMasterFacesPerEdge, edgeI)
4179 if (nMasterFacesPerEdge[edgeI] > 2)
4186 Info<<
"Protected from visiting " 4188 <<
" non-manifold edges" <<
nl <<
endl;
4193 DynamicList<label> changedEdges;
4194 DynamicList<patchFaceOrientation> changedInfo;
4196 const scalar tol = PatchEdgeFaceWave
4199 patchFaceOrientation
4200 >::propagationTol();
4204 globalIndex globalFaces(
patch.size());
4214 globalSeed = globalFaces.toGlobal(faceI);
4219 reduce(globalSeed, minOp<label>());
4226 if (globalFaces.isLocal(globalSeed))
4228 const label seedFaceI = globalFaces.toLocal(globalSeed);
4232 patchFaceOrientation& faceInfo =
allFaceInfo[seedFaceI];
4237 if (zoneToOrientation[faceToZone[seedFaceI]] < 0)
4246 label edgeI = fEdges[fEdgeI];
4248 patchFaceOrientation& edgeInfo = allEdgeInfo[edgeI];
4252 edgeInfo.updateEdge<
int>
4264 changedEdges.append(edgeI);
4265 changedInfo.append(edgeInfo);
4282 patchFaceOrientation
4300 mesh_.nBoundaryFaces(),
4306 const label meshFaceI =
patch.addressing()[i];
4307 if (!mesh_.isInternalFace(meshFaceI))
4309 neiStatus[meshFaceI-mesh_.nInternalFaces()] =
4317 const label meshFaceI =
patch.addressing()[i];
4318 const label patchI = bm.whichPatch(meshFaceI);
4324 && !isMasterFace.test(meshFaceI)
4328 label bFaceI = meshFaceI-mesh_.nInternalFaces();
4341 <<
"Incorrect status for face " << meshFaceI
4351 meshFlipMap.setSize(mesh_.nFaces());
4352 meshFlipMap =
false;
4356 label meshFaceI =
patch.addressing()[faceI];
4360 meshFlipMap.unset(meshFaceI);
4364 meshFlipMap.set(meshFaceI);
4369 <<
"Problem : unvisited face " << faceI
4370 <<
" centre:" << mesh_.faceCentres()[meshFaceI]
4377 void Foam::meshRefinement::zonify
4380 const bitSet& isMasterFace,
4384 const bitSet& meshFlipMap,
4385 polyTopoChange& meshMod
4388 const labelList& faceOwner = mesh_.faceOwner();
4389 const labelList& faceNeighbour = mesh_.faceNeighbour();
4391 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
4393 label faceZoneI = faceToZone[faceI];
4395 if (faceZoneI != -1)
4401 label ownZone = cellToZone[faceOwner[faceI]];
4402 label neiZone = cellToZone[faceNeighbour[faceI]];
4406 if (ownZone == neiZone)
4409 flip = meshFlipMap.
test(faceI);
4416 || (neiZone != -1 && ownZone > neiZone)
4424 mesh_.faces()[faceI],
4427 faceNeighbour[faceI],
4439 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
4446 label faceI =
pp.start();
4450 label faceZoneI = faceToZone[faceI];
4452 if (faceZoneI != -1)
4454 label ownZone = cellToZone[faceOwner[faceI]];
4455 label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
4459 if (ownZone == neiZone)
4462 flip = meshFlipMap.test(faceI);
4469 || (neiZone != -1 && ownZone > neiZone)
4477 mesh_.faces()[faceI],
4497 forAll(cellToZone, cellI)
4499 label zoneI = cellToZone[cellI];
4517 void Foam::meshRefinement::allocateInterRegionFaceZone
4519 const label ownZone,
4520 const label neiZone,
4522 LabelPairMap<word>& zoneIDsToFaceZone
4527 if (ownZone != neiZone)
4534 || (neiZone != -1 && ownZone > neiZone)
4544 if (!zoneIDsToFaceZone.found(
key))
4547 const word ownZoneName =
4550 ? cellZones[ownZone].name()
4553 const word neiZoneName =
4556 ? cellZones[neiZone].name()
4561 Pair<word> wordKey(ownZoneName, neiZoneName);
4567 word fzName = wordKey.first() +
"_to_" + wordKey.second();
4569 zoneIDsToFaceZone.insert(
key, fzName);
4570 zonesToFaceZone.insert(wordKey, fzName);
4580 const bool doHandleSnapProblems,
4582 const bool useTopologicalSnapDetection,
4583 const bool removeEdgeConnectedCells,
4585 const label nErodeCellZones,
4594 const bool exitIfLeakPath,
4605 Info<<
"Introducing baffles for " 4607 <<
" faces that are intersected by the surface." <<
nl <<
endl;
4610 labelList neiLevel(mesh_.nBoundaryFaces());
4612 calcNeighbourData(neiLevel, neiCc);
4620 globalToMasterPatch,
4624 locationsOutsideMesh,
4636 createBaffles(ownPatch, neiPatch);
4644 Info<<
"Created baffles in = " 4647 printMeshInfo(
debug,
"After introducing baffles",
true);
4651 const_cast<Time&
>(mesh_.time())++;
4660 Pout<<
"Dumped debug data in = " 4670 if (doHandleSnapProblems)
4675 useTopologicalSnapDetection,
4676 removeEdgeConnectedCells,
4680 globalToMasterPatch,
4688 neiLevel.setSize(mesh_.nBoundaryFaces());
4689 neiCc.setSize(mesh_.nBoundaryFaces());
4690 calcNeighbourData(neiLevel, neiCc);
4697 globalToMasterPatch,
4701 locationsOutsideMesh,
4713 createBaffles(ownPatch, neiPatch);
4728 <<
"Remove unreachable sections of mesh" <<
nl 4729 <<
"-----------------------------------" <<
nl 4739 globalToMasterPatch,
4742 locationsOutsideMesh,
4752 Info<<
"Split mesh in = " 4755 printMeshInfo(
debug,
"After subsetting",
true);
4768 Pout<<
"Dumped debug data in = " 4777 const bool useTopologicalSnapDetection,
4778 const bool removeEdgeConnectedCells,
4780 const scalar planarAngle,
4793 <<
"Merge free-standing baffles" <<
nl 4794 <<
"---------------------------" <<
nl 4810 Info<<
"Detected free-standing baffles : " << nCouples <<
endl;
4824 useTopologicalSnapDetection,
4825 removeEdgeConnectedCells,
4829 globalToMasterPatch,
4837 <<
"Remove unreachable sections of mesh" <<
nl 4838 <<
"-----------------------------------" <<
nl 4848 globalToMasterPatch,
4851 locationsOutsideMesh,
4863 Info<<
"Merged free-standing baffles in = " 4870 const label nBufferLayers,
4871 const label nErodeCellZones,
4878 const bool exitIfLeakPath,
4888 labelList neiLevel(mesh_.nBoundaryFaces());
4890 calcNeighbourData(neiLevel, neiCc);
4897 globalToMasterPatch,
4901 locationsOutsideMesh,
4914 boolList blockedFace(mesh_.nFaces(),
false);
4918 if (ownPatch[faceI] != -1 || neiPatch[faceI] != -1)
4920 blockedFace[faceI] =
true;
4932 vector::uniform(mergeDistance_),
4934 locationsOutsideMesh,
4946 globalToMasterPatch,
4957 const label nBufferLayers,
4972 const labelList& faceOwner = mesh_.faceOwner();
4973 const labelList& faceNeighbour = mesh_.faceNeighbour();
4977 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
4979 if (ownPatch[facei] == -1 && neiPatch[facei] != -1)
4982 <<
" at:" << mesh_.faceCentres()[facei]
4983 <<
" ownPatch:" << ownPatch[facei]
4984 <<
" neiPatch:" << neiPatch[facei]
4990 const label ownRegion = cellRegion[faceOwner[facei]];
4991 const label neiRegion = cellRegion[faceNeighbour[facei]];
4992 if (ownRegion != neiRegion)
4994 if (ownPatch[facei] == -1)
4997 <<
" at:" << mesh_.faceCentres()[facei]
4998 <<
" ownPatch:" << ownPatch[facei]
4999 <<
" ownRegion:" << ownRegion
5000 <<
" neiPatch:" << neiPatch[facei]
5001 <<
" neiRegion:" << neiRegion
5013 DynamicList<label> startFaces;
5016 if (ownPatch[facei] != -1)
5018 startFaces.append(facei);
5024 autoPtr<mapDistribute> mapPtr;
5028 bitSet(mesh_.nFaces()),
5036 labelList startOwnPatch(ownPatch, startFaces);
5037 mapPtr().distribute(startOwnPatch);
5039 nearestOwnPatch.setSize(mesh_.nFaces());
5040 nearestOwnPatch = -1;
5041 forAll(faceToStart, facei)
5043 const label sloti = faceToStart[facei];
5046 nearestOwnPatch[facei] = startOwnPatch[sloti];
5060 bitSet isFrozenPoint(mesh_.nPoints());
5061 bitSet isFrozenFace(mesh_.nFaces());
5067 const label pointZonei = pzs.
findZoneID(
"frozenPoints");
5068 if (pointZonei != -1)
5070 const pointZone& pz = pzs[pointZonei];
5071 isFrozenPoint.set(pz);
5072 for (
const label pointi : pz)
5074 isFrozenFace.set(pointFaces[pointi]);
5079 const label faceZonei = fzs.
findZoneID(
"frozenFaces");
5080 if (faceZonei != -1)
5082 const faceZone& fz = fzs[faceZonei];
5083 isFrozenFace.set(fz);
5084 for (
const label facei : fz)
5086 isFrozenPoint.set(mesh_.faces()[facei]);
5092 label defaultPatch = 0;
5093 if (globalToMasterPatch.
size())
5095 defaultPatch = globalToMasterPatch[0];
5098 for (label i = 0; i < nBufferLayers; i++)
5102 labelList pointBaffle(mesh_.nPoints(), -1);
5104 forAll(faceNeighbour, facei)
5106 if (!isFrozenFace[facei])
5108 const face&
f = mesh_.faces()[facei];
5110 const label ownRegion = cellRegion[faceOwner[facei]];
5111 const label neiRegion = cellRegion[faceNeighbour[facei]];
5113 if (ownRegion == -1 && neiRegion != -1)
5120 if (!isFrozenPoint[
f[fp]])
5122 pointBaffle[
f[fp]] =
5123 max(defaultPatch, ownPatch[facei]);
5127 else if (ownRegion != -1 && neiRegion == -1)
5129 label newPatchi = neiPatch[facei];
5130 if (newPatchi == -1)
5132 newPatchi =
max(defaultPatch, ownPatch[facei]);
5136 if (!isFrozenPoint[
f[fp]])
5138 pointBaffle[
f[fp]] = newPatchi;
5177 const auto&
patches = mesh_.boundaryMesh();
5189 const label facei =
pp.start()+i;
5190 if (!isFrozenFace[facei])
5192 const face&
f = mesh_.faces()[facei];
5193 const label ownRegion = cellRegion[faceOwner[facei]];
5194 const label neiRegion =
5195 neiCellRegion[facei-mesh_.nInternalFaces()];
5199 if (ownRegion == -1 && neiRegion != -1)
5203 if (!isFrozenPoint[
f[fp]])
5205 pointBaffle[
f[fp]] =
5206 max(defaultPatch, ownPatch[facei]);
5210 else if (ownRegion != -1 && neiRegion == -1)
5212 label newPatchI = neiPatch[facei];
5213 if (newPatchI == -1)
5215 newPatchI =
max(defaultPatch, ownPatch[facei]);
5219 if (!isFrozenPoint[
f[fp]])
5221 pointBaffle[
f[fp]] = newPatchI;
5232 const label facei =
pp.start()+i;
5233 if (!isFrozenFace[facei])
5235 const face&
f = mesh_.faces()[facei];
5236 const label ownRegion = cellRegion[faceOwner[facei]];
5238 if (ownRegion != -1)
5242 if (!isFrozenPoint[
f[fp]])
5244 pointBaffle[
f[fp]] =
5245 max(defaultPatch, ownPatch[facei]);
5269 forAll(pointFaces, pointi)
5271 if (pointBaffle[pointi] != -1)
5277 const label facei =
pFaces[pFacei];
5279 if (!isFrozenFace[facei] && ownPatch[facei] == -1)
5281 ownPatch[facei] = pointBaffle[pointi];
5295 if (!isFrozenFace[facei] && ownPatch[facei] != -1)
5297 const label own = faceOwner[facei];
5299 if (cellRegion[own] == -1)
5303 const cell& ownFaces = mesh_.cells()[own];
5306 const label ownFacei = ownFaces[j];
5307 if (!isFrozenFace[ownFacei] && ownPatch[ownFacei] == -1)
5309 newOwnPatch[ownFacei] = ownPatch[facei];
5313 if (mesh_.isInternalFace(facei))
5315 const label nei = faceNeighbour[facei];
5317 if (cellRegion[nei] == -1)
5321 const cell& neiFaces = mesh_.cells()[nei];
5324 const label neiFacei = neiFaces[j];
5325 const bool isFrozen = isFrozenFace[neiFacei];
5326 if (!isFrozen && ownPatch[neiFacei] == -1)
5328 newOwnPatch[neiFacei] = ownPatch[facei];
5346 DynamicList<label> cellsToRemove(mesh_.nCells());
5347 forAll(cellRegion, celli)
5349 if (cellRegion[celli] == -1)
5351 cellsToRemove.append(celli);
5354 cellsToRemove.shrink();
5358 mesh_.nCells() - cellsToRemove.size(),
5362 Info<<
"Keeping all cells containing inside points" <<
endl 5363 <<
"Selected for keeping : " << nCellsToKeep <<
" cells." <<
endl;
5367 removeCells cellRemover(mesh_);
5370 const labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
5371 labelList exposedPatches(exposedFaces.size());
5373 label nUnpatched = 0;
5377 label facei = exposedFaces[i];
5379 if (ownPatch[facei] != -1)
5381 exposedPatches[i] = ownPatch[facei];
5385 const label fallbackPatch =
5387 nearestOwnPatch.size()
5388 ? nearestOwnPatch[facei]
5391 if (nUnpatched == 0)
5394 <<
"For exposed face " << facei
5395 <<
" fc:" << mesh_.faceCentres()[facei]
5396 <<
" found no patch." <<
endl 5397 <<
" Taking patch " << fallbackPatch
5398 <<
" instead. Suppressing future warnings" <<
endl;
5402 exposedPatches[i] = fallbackPatch;
5406 reduce(nUnpatched, sumOp<label>());
5409 Info<<
"Detected " << nUnpatched <<
" faces out of " 5411 <<
" for which the default patch " << defaultPatch
5412 <<
" will be used" <<
endl;
5415 return doRemoveCells
5427 const label nBufferLayers,
5428 const label nErodeCellZones,
5440 labelList neiLevel(mesh_.nBoundaryFaces());
5442 calcNeighbourData(neiLevel, neiCc);
5450 globalToMasterPatch,
5454 locationsOutsideMesh,
5470 limitShells_.findLevel
5472 mesh_.cellCentres(),
5476 forAll(levelShell, celli)
5478 if (levelShell[celli] != -1)
5481 cellRegion[celli] = -1;
5488 globalToMasterPatch,
5497 const_cast<Time&
>(mesh_.time())++;
5498 Pout<<
"Writing mesh after removing limitShells" 5529 Info<<
"dupNonManifoldPoints : Found : " << nNonManifPoints
5530 <<
" non-manifold points (out of " 5531 << mesh_.globalData().nTotalPoints()
5537 if (nNonManifPoints)
5547 mesh_.moving(
false);
5550 mapPtr = meshMod.
changeMesh(mesh_,
false,
true);
5554 mesh_.updateMesh(map);
5593 label nPointPairs = 0;
5594 forAll(pointToDuplicate, pointI)
5596 label otherPointI = pointToDuplicate[pointI];
5597 if (otherPointI != -1)
5608 forAll(pointToDuplicate, pointI)
5610 label otherPointI = pointToDuplicate[pointI];
5611 if (otherPointI != -1)
5614 pointToMaster.
insert(pointI, otherPointI);
5626 mesh_.moving(
false);
5629 mapPtr = meshMod.
changeMesh(mesh_,
false,
true);
5633 mesh_.updateMesh(map);
5680 internalOrBaffleFaceZones = getZones(fzTypes);
5690 forAll(boundaryFaceZones, j)
5692 const faceZone& fZone = mesh_.faceZones()[boundaryFaceZones[j]];
5695 const face&
f = mesh_.faces()[fZone[i]];
5698 pointStatus[
f[fp]] =
max(pointStatus[
f[fp]], 1u);
5702 forAll(internalOrBaffleFaceZones, j)
5704 const faceZone& fZone = mesh_.faceZones()[internalOrBaffleFaceZones[j]];
5707 const face&
f = mesh_.faces()[fZone[i]];
5710 pointStatus[
f[fp]] =
max(pointStatus[
f[fp]], 2u);
5725 forAll(pointStatus, pointI)
5727 if (pointStatus[pointI] == 1u)
5734 Info<<
"Duplicating " << globalNPoints <<
" points on" 5735 <<
" faceZones of type " 5745 forAll(pointStatus, pointI)
5747 if (pointStatus[pointI] == 1u)
5749 candidatePoints[
n++] = pointI;
5762 const bool allowFreeStandingZoneFaces,
5763 const label nErodeCellZones,
5767 const bool exitIfLeakPath,
5773 if (locationsInMesh.
size() != zonesInMesh.
size())
5783 labelList neiLevel(mesh_.nBoundaryFaces());
5785 calcNeighbourData(neiLevel, neiCc);
5794 if (namedSurfaces.
size())
5796 Info<<
"Setting cellZones according to named surfaces:" <<
endl;
5799 label surfI = namedSurfaces[i];
5801 Info<<
"Surface : " << surfaces_.names()[surfI] <<
nl 5802 <<
" faceZones : " << surfZones[surfI].faceZoneNames() <<
nl 5803 <<
" cellZone : " << surfZones[surfI].cellZoneName()
5844 allowFreeStandingZoneFaces,
5849 locationsOutsideMesh,
5867 labelList faceToZone(mesh_.nFaces(), -1);
5869 forAll(namedSurfaceRegion, faceI)
5872 label globalI = namedSurfaceRegion[faceI];
5875 const labelPair spr = surfaces_.whichSurface(globalI);
5876 faceToZone[faceI] = surfaceToFaceZones[spr.
first()][spr.
second()];
5890 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
5892 if (faceToZone[faceI] == -1)
5897 allocateInterRegionFaceZone
5899 cellToZone[mesh_.faceOwner()[faceI]],
5900 cellToZone[mesh_.faceNeighbour()[faceI]],
5910 forAll(neiCellZone, bFaceI)
5912 label faceI = bFaceI + mesh_.nInternalFaces();
5913 if (faceToZone[faceI] == -1)
5915 allocateInterRegionFaceZone
5917 cellToZone[mesh_.faceOwner()[faceI]],
5918 neiCellZone[bFaceI],
5936 Info<<
"Setting faceZones according to neighbouring cellZones:" 5950 const word& fzName = zonesToFaceZone[cz];
5953 << cz[0] <<
' ' << cz[1] <<
nl 5954 <<
" faceZone : " << fzName <<
endl;
5964 label cz0 = cellZones.findZoneID(cz[0]);
5965 label cz1 = cellZones.findZoneID(cz[1]);
5975 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
5977 if (faceToZone[faceI] == -1)
5983 label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
5984 label neiZone = cellToZone[mesh_.faceNeighbour()[faceI]];
5985 if (ownZone != neiZone)
5990 || (neiZone != -1 && ownZone > neiZone)
5997 faceToZone[faceI] = fZoneLookup[
key];
6001 forAll(neiCellZone, bFaceI)
6003 label faceI = bFaceI + mesh_.nInternalFaces();
6004 if (faceToZone[faceI] == -1)
6006 label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
6007 label neiZone = neiCellZone[bFaceI];
6008 if (ownZone != neiZone)
6013 || (neiZone != -1 && ownZone > neiZone)
6020 faceToZone[faceI] = fZoneLookup[
key];
6039 label bFaceI =
pp.start()-mesh_.nInternalFaces();
6042 neiCellZone[bFaceI++] = -1;
6063 bitSet meshFlipMap(mesh_.nFaces(),
false);
6071 freeStandingBaffleFaces
6082 if (nFreeStanding > 0)
6084 Info<<
"Detected " << nFreeStanding <<
" free-standing zone faces" 6089 OBJstream str(mesh_.time().path()/
"freeStanding.obj");
6090 Pout<<
"meshRefinement::zonify : dumping faceZone faces to " 6098 calcPatchNumMasterFaces(isMasterFace,
patch, nMasterFacesPerEdge);
6103 const label nZones = markPatchZones
6106 nMasterFacesPerEdge,
6111 for (label zoneI = 0; zoneI < nZones; zoneI++)
6113 nPosOrientation.
insert(zoneI, 0);
6119 consistentOrientation
6123 nMasterFacesPerEdge,
6124 faceToConnectedZone,
6134 label meshFaceI =
patch.addressing()[faceI];
6136 if (isMasterFace.
test(meshFaceI))
6141 posOrientation.
test(meshFaceI)
6142 == meshFlipMap.test(meshFaceI)
6148 nPosOrientation.
find(faceToConnectedZone[faceI])() +=
n;
6154 Info<<
"Split " << nFreeStanding <<
" free-standing zone faces" 6155 <<
" into " << nZones <<
" disconnected regions with size" 6156 <<
" (negative denotes wrong orientation) :" 6159 for (label zoneI = 0; zoneI < nZones; zoneI++)
6161 Info<<
" " << zoneI <<
"\t" << nPosOrientation[zoneI]
6169 consistentOrientation
6173 nMasterFacesPerEdge,
6174 faceToConnectedZone,
6201 mesh_.moving(
false);
6207 mesh_.updateMesh(map());
6210 if (map().hasMotionPoints())
6212 mesh_.movePoints(map().preMotionPoints());
6224 if (mesh_.cellZones().size() > 0)
6227 forAll(mesh_.cellZones(), zoneI)
6229 const cellZone& cz = mesh_.cellZones()[zoneI];
6236 if (mesh_.faceZones().size() > 0)
6239 forAll(mesh_.faceZones(), zoneI)
6241 const faceZone& fz = mesh_.faceZones()[zoneI];
static const word & zeroGradientType() noexcept
The type name for zeroGradient patch fields.
This class separates the mesh into distinct unconnected regions, each of which is then given a label ...
Class containing data for face removal.
const polyBoundaryMesh & pbm
const T & first() const noexcept
Access the first element.
void size(const label n)
Older name for setAddressableSize.
ZoneMesh< faceZone, polyMesh > faceZoneMesh
A ZoneMesh with the type faceZone.
const labelIOList & zoneIDs
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
fileName path() const
Return path = rootPath/caseName. Same as TimePaths::path()
labelList findIndices(const ListType &input, const UnaryPredicate &pred, label start=0)
Linear search to find all occurences of given element.
Ostream & indent(Ostream &os)
Indent stream.
List< wallPoints > allCellInfo(mesh_.nCells())
autoPtr< mapPolyMesh > mergePoints(const labelList &pointToDuplicate)
Merge duplicate points.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Class describing modification of a face.
A face is a list of labels corresponding to mesh vertices.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void transfer(List< T > &list)
Transfer the contents of the argument List into this list and annul the argument list.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
void append(const T &val)
Append an element at the end of the list.
virtual const labelList & faceNeighbour() const
Return face neighbour.
const bitSet isBlockedFace(intersectedFaces())
wallPoints::trackData td(isBlockedFace, regionToBlockSize)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Unit conversion functions.
constexpr char nl
The newline '\n' character (0x0a)
static labelListList addFaceZonesToMesh(const PtrList< surfaceZonesInfo > &surfList, const labelList &namedSurfaces, polyMesh &mesh)
UIndirectList< label > labelUIndList
UIndirectList of labels.
T & first()
Access first element of the list, position [0].
Ostream & endl(Ostream &os)
Add newline and flush stream.
bool found(const T &val, label pos=0) const
Same as contains()
autoPtr< mapPolyMesh > createBaffles(const labelList &ownPatch, const labelList &neiPatch)
Create baffle for every internal face where ownPatch != -1.
void baffleAndSplitMesh(const bool handleSnapProblems, const snapParameters &snapParams, const bool useTopologicalSnapDetection, const bool removeEdgeConnectedCells, const scalarField &perpendicularAngle, const label nErodeCellZones, const dictionary &motionDict, Time &runTime, const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, const pointField &locationsInMesh, const wordList ®ionsInMesh, const pointField &locationsOutsideMesh, const bool exitIfLeakPath, const refPtr< coordSetWriter > &leakPathFormatter, refPtr< surfaceWriter > &surfFormatter)
Split off unreachable areas of mesh.
static labelList addCellZonesToMesh(const PtrList< surfaceZonesInfo > &surfList, const labelList &namedSurfaces, polyMesh &mesh)
static const Enum< faceZoneType > faceZoneTypeNames
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
labelList getZones(const List< surfaceZonesInfo::faceZoneType > &fzTypes) const
Get zones of given type.
static List< labelPair > subsetBaffles(const polyMesh &mesh, const labelList &zoneIDs, const List< labelPair > &baffles)
Subset baffles according to zones.
Ignore writing from objectRegistry::writeObject()
const dimensionSet dimless
Dimensionless.
label nFaces() const noexcept
Number of mesh faces.
static labelList getStandaloneNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of named surfaces without a cellZone.
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
A PrimitivePatch with an IndirectList for the faces, const reference for the point field...
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.
virtual const fileName & name() const override
Get the name of the output serial stream. (eg, the name of the Fstream file name) ...
labelList faceLabels(nFaceLabels)
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
List< labelList > labelListList
List of labelList.
DynamicID< faceZoneMesh > faceZoneID
Foam::faceZoneID.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
A class for managing references or pointers (no reference counting)
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
Takes mesh with 'baffles' (= boundary faces sharing points). Determines for selected points on bounda...
static labelPairList findDuplicateFacePairs(const polyMesh &)
Helper routine to find all baffles (two boundary faces.
#define forAll(list, i)
Loop across all elements in list.
const labelList & reverseFaceMap() const noexcept
Reverse face map.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
bool insert(const Key &key, const T &obj)
Copy insert a new entry, not overwriting existing entries.
void mergeFreeStandingBaffles(const snapParameters &snapParams, const bool useTopologicalSnapDetection, const bool removeEdgeConnectedCells, const scalarField &perpendicularAngle, const scalar planarAngle, const dictionary &motionDict, Time &runTime, const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, const pointField &locationsInMesh, const pointField &locationsOutsideMesh)
Merge free-standing baffles.
label size() const noexcept
The number of elements in table.
void write(vtk::formatter &fmt, const Type &val, const label n=1)
Component-wise write of a value (N times)
A list of faces which address into the list of points.
Determines the 'side' for every face and connected to a singly-connected (through edges) region of fa...
bool returnReduceAnd(const bool value, const label comm=UPstream::worldComm)
Perform logical (and) MPI Allreduce on a copy. Uses UPstream::reduceAnd.
word outputName("finiteArea-edges.obj")
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
void setSize(const label n)
Alias for resize()
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
dimensionedScalar cos(const dimensionedScalar &ds)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
label findZoneID(const word &zoneName) const
Find zone index by name, return -1 if not found.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i), works like std::iota() but returning a...
static void mergePoints(const polyMesh &, const Map< label > &pointToMaster, polyTopoChange &meshMod)
Helper: Merge points.
const wordList area
Standard area field types (scalar, vector, tensor, etc)
void clear()
Clear the list, i.e. set size to zero.
A class for handling words, derived from Foam::string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
iterator find(const Key &key)
Find and return an iterator set at the hashed entry.
Simple container to keep together snap specific information.
label size() const noexcept
The number of entries in the list.
const Field< point_type > & faceCentres() const
Return face centres for patch.
autoPtr< mapPolyMesh > mergeZoneBaffles(const bool doInternalZones, const bool doBaffleZones)
Merge all baffles on faceZones.
virtual const labelList & faceOwner() const
Return face owner.
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
autoPtr< mapPolyMesh > removeLimitShells(const label nBufferLayers, const label nErodeCellZones, const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, const pointField &locationsInMesh, const wordList ®ionsInMesh, const pointField &locationsOutsideMesh)
Remove cells from limitRegions if level -1.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
ZoneMesh< pointZone, polyMesh > pointZoneMesh
A ZoneMesh with the type pointZone.
void append(const T &val)
Copy append an element to the end of this list.
virtual const faceList & faces() const
Return raw faces.
label whichPatch(const label meshFacei) const
Return patch index for a given mesh face index. Uses binary search.
static labelList getInsidePointNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of surfaces with a cellZone that have 'insidePoint'.
FaceCellWave< wallPoints, wallPoints::trackData > wallDistCalc(mesh_, changedFaces, faceDist, allFaceInfo, allCellInfo, 0, td)
A HashTable similar to std::unordered_map.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
static labelList getUnnamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of unnamed surfaces (surfaces without faceZoneName)
errorManip< error > abort(error &err)
const T * set(const label i) const
Return const pointer to element (can be nullptr), or nullptr for out-of-range access (ie...
bool test(const label pos) const
Test for True value at specified position, never auto-vivify entries.
label find(const T &val) const
Find index of the first occurrence of the value.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO...
Istream and Ostream manipulators taking arguments.
An OFstream that keeps track of vertices and provides convenience output methods for OBJ files...
Smanip< std::ios_base::fmtflags > setf(std::ios_base::fmtflags flags)
int debug
Static debugging option.
Pair< label > labelPair
A pair of labels.
Type gMax(const FieldField< Field, Type > &f)
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
const faceZoneMesh & faceZones() const noexcept
Return face zone mesh.
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type test(const label i) const
Test bool value at specified position, always false for out-of-range access.
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]
double cpuTimeIncrement() const
Return CPU time [seconds] since last call to cpuTimeIncrement(), resetCpuTimeIncrement().
bool hasMotionPoints() const noexcept
Has valid preMotionPoints?
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
autoPtr< mapPolyMesh > createZoneBaffles(const labelList &zoneIDs, List< labelPair > &baffles, labelList &originatingFaceZone)
Create baffles for faces on faceZones. Return created baffles.
const vectorField & faceCentres() const
List< word > wordList
List of word.
PrimitivePatch< UIndirectList< face >, const pointField & > uindirectPrimitivePatch
A PrimitivePatch with UIndirectList for the faces, const reference for the point field.
autoPtr< mapPolyMesh > dupNonManifoldPoints()
Find boundary points that connect to more than one cell.
vector point
Point is a vector.
autoPtr< mapPolyMesh > dupNonManifoldBoundaryPoints()
Find boundary points that are on faceZones of type boundary.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
#define WarningInFunction
Report a warning using Foam::Warning.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
auto key(const Type &t) -> typename std::enable_if< std::is_enum< Type >::value, typename std::underlying_type< Type >::type >::type
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers...
List< Key > sortedToc() const
The table of contents (the keys) in sorted order.
Direct mesh changes based on v1.3 polyTopoChange syntax.
const polyBoundaryMesh & patches
const word & name() const noexcept
The zone name.
void getZoneFaces(const labelList &zoneIDs, labelList &faceZoneID, labelList &ownPatch, labelList &neiPatch, labelList &nBaffles) const
Get per-face information (faceZone, master/slave patch)
static word outputPrefix
Directory prefix.
const std::string patch
OpenFOAM patch number as a std::string.
void setRefinement(const localPointRegion ®ionSide, polyTopoChange &)
Play commands into polyTopoChange to duplicate points. Gets.
label nRegions() const
Return total number of regions.
static List< pointField > zonePoints(const pointField &locationsInMesh, const wordList &zonesInMesh, const pointField &locationsOutsideMesh)
Helper: per zone (entry in zonesInMesh) the locations with.
autoPtr< mapPolyMesh > mergeBaffles(const List< labelPair > &, const Map< label > &faceToPatch)
Merge baffles. Gets pairs of faces and boundary faces to move.
static labelList getClosedNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList, const searchableSurfaces &allGeometry, const labelList &surfaces)
Get indices of surfaces with a cellZone that are closed and.
const boolList & flipMap() const noexcept
Return face flip map.
static labelList getNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of named surfaces (surfaces with faceZoneName)
messageStream Info
Information stream (stdout output on master, null elsewhere)
writeType
Enumeration for what to write. Used as a bit-pattern.
Field< vector > vectorField
Specialisation of Field<T> for vector.
static autoPtr< mapDistribute > calcClosure(const polyMesh &mesh, const List< pointField > &zonePoints, const labelList &blockedFaces, const globalIndex &globalBlockedFaces, const bool erode, labelList &closureFaces, labelList &closureToBlocked)
Optional direct use to generate the set of faces and the method to.
faceZoneType
What to do with faceZone faces.
debugType
Enumeration for what to debug. Used as a bit-pattern.
Mesh consisting of general polyhedral cells.
Omanip< int > setw(const int i)
A subset of mesh faces organised as a primitive patch.
List< label > labelList
A List of labels.
A patch is a list of labels that address the faces in the global face list.
label checkData(const objectRegistry &obr, const instantList &timeDirs, wordList &objectNames, const fileName &local=fileName::null)
Check if fields are good to use (available at all times)
const T & second() const noexcept
Access the second element.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
static void mapCombineReduce(Container &values, const CombineOp &cop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce) applying cop to inplace combine map values from different processo...
bool returnReduceOr(const bool value, const label comm=UPstream::worldComm)
Perform logical (or) MPI Allreduce on a copy. Uses UPstream::reduceOr.
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
List< bool > boolList
A List of bools.
A List with indirect addressing.
const pointField & preMotionPoints() const noexcept
Pre-motion point positions.
List< wallPoints > allFaceInfo(mesh_.nFaces())
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Do not request registration (bool: false)
HashTable< word, wordPair, Foam::Hash< wordPair > > wordPairHashTable
HashTable of wordPair.
ZoneMesh< cellZone, polyMesh > cellZoneMesh
A ZoneMesh with the type cellZone.
void reduce(T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce) using linear/tree communication schedule.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
static label addFaceZone(const word &name, const labelList &addressing, const boolList &flipMap, polyMesh &mesh)
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
forAllConstIters(mixture.phases(), phase)
static void listCombineReduce(List< T > &values, const CombineOp &cop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Combines List elements. After completion all processors have the same data.
A HashTable to objects of type <T> with a label key.
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.
static constexpr const zero Zero
Global zero (0)