134 using namespace Foam;
138 template<
class GeoField>
143 mesh.objectRegistry::lookupClass<GeoField>()
153 const GeoField& fld = *iter();
155 typename GeoField::GeometricBoundaryField& bfld =
156 const_cast<typename GeoField::GeometricBoundaryField&
>
161 label sz = bfld.size();
166 GeoField::PatchFieldType::New
170 fld.dimensionedInternalField()
178 template<
class GeoField>
183 mesh.objectRegistry::lookupClass<GeoField>()
193 const GeoField& fld = *iter();
195 const_cast<typename GeoField::GeometricBoundaryField&
>
204 template<
class GeoField>
209 mesh.objectRegistry::lookupClass<GeoField>()
219 const GeoField& fld = *iter();
221 typename GeoField::GeometricBoundaryField& bfld =
222 const_cast<typename GeoField::GeometricBoundaryField&
>
227 bfld.reorder(oldToNew);
241 if (polyPatches[patchI].
type() == patch.type())
249 <<
"Already have patch " << patch.
name()
250 <<
" but of type " << patch.type()
256 label insertPatchI = polyPatches.
size();
257 label startFaceI = mesh.
nFaces();
259 forAll(polyPatches, patchI)
261 const polyPatch& pp = polyPatches[patchI];
263 if (isA<processorPolyPatch>(pp))
265 insertPatchI = patchI;
266 startFaceI = pp.
start();
278 label sz = polyPatches.
size();
306 addPatchFields<volScalarField>
311 addPatchFields<volVectorField>
316 addPatchFields<volSphericalTensorField>
321 addPatchFields<volSymmTensorField>
326 addPatchFields<volTensorField>
334 addPatchFields<surfaceScalarField>
339 addPatchFields<surfaceVectorField>
344 addPatchFields<surfaceSphericalTensorField>
349 addPatchFields<surfaceSymmTensorField>
354 addPatchFields<surfaceTensorField>
363 for (label i = 0; i < insertPatchI; i++)
368 for (label i = insertPatchI; i < sz; i++)
373 oldToNew[sz] = insertPatchI;
379 reorderPatchFields<volScalarField>(
mesh, oldToNew);
380 reorderPatchFields<volVectorField>(
mesh, oldToNew);
381 reorderPatchFields<volSphericalTensorField>(
mesh, oldToNew);
382 reorderPatchFields<volSymmTensorField>(
mesh, oldToNew);
383 reorderPatchFields<volTensorField>(
mesh, oldToNew);
384 reorderPatchFields<surfaceScalarField>(
mesh, oldToNew);
385 reorderPatchFields<surfaceVectorField>(
mesh, oldToNew);
386 reorderPatchFields<surfaceSphericalTensorField>(
mesh, oldToNew);
387 reorderPatchFields<surfaceSymmTensorField>(
mesh, oldToNew);
388 reorderPatchFields<surfaceTensorField>(
mesh, oldToNew);
399 const label nNewPatches
410 reorderPatchFields<volScalarField>(
mesh, oldToNew);
411 reorderPatchFields<volVectorField>(
mesh, oldToNew);
412 reorderPatchFields<volSphericalTensorField>(
mesh, oldToNew);
413 reorderPatchFields<volSymmTensorField>(
mesh, oldToNew);
414 reorderPatchFields<volTensorField>(
mesh, oldToNew);
415 reorderPatchFields<surfaceScalarField>(
mesh, oldToNew);
416 reorderPatchFields<surfaceVectorField>(
mesh, oldToNew);
417 reorderPatchFields<surfaceSphericalTensorField>(
mesh, oldToNew);
418 reorderPatchFields<surfaceSymmTensorField>(
mesh, oldToNew);
419 reorderPatchFields<surfaceTensorField>(
mesh, oldToNew);
422 polyPatches.
setSize(nNewPatches);
423 fvPatches.
setSize(nNewPatches);
424 trimPatchFields<volScalarField>(
mesh, nNewPatches);
425 trimPatchFields<volVectorField>(
mesh, nNewPatches);
426 trimPatchFields<volSphericalTensorField>(
mesh, nNewPatches);
427 trimPatchFields<volSymmTensorField>(
mesh, nNewPatches);
428 trimPatchFields<volTensorField>(
mesh, nNewPatches);
429 trimPatchFields<surfaceScalarField>(
mesh, nNewPatches);
430 trimPatchFields<surfaceVectorField>(
mesh, nNewPatches);
431 trimPatchFields<surfaceSphericalTensorField>(
mesh, nNewPatches);
432 trimPatchFields<surfaceSymmTensorField>(
mesh, nNewPatches);
433 trimPatchFields<surfaceTensorField>(
mesh, nNewPatches);
437 template<
class GeoField>
451 mesh.objectRegistry::lookupClass<GeoField>()
455 const GeoField& fld = *iter();
457 Info<<
"Mapping field " << fld.name() <<
endl;
475 if (addedPatches.
found(patchI))
477 tSubFld().boundaryField()[patchI] ==
483 GeoField* subFld = tSubFld.ptr();
484 subFld->rename(fld.name());
491 template<
class GeoField>
492 void subsetSurfaceFields
504 mesh.objectRegistry::lookupClass<GeoField>()
508 const GeoField& fld = *iter();
510 Info<<
"Mapping field " << fld.name() <<
endl;
527 if (addedPatches.
found(patchI))
529 tSubFld().boundaryField()[patchI] ==
535 GeoField* subFld = tSubFld.ptr();
536 subFld->rename(fld.name());
548 if (cellRegion[cellI] != regionI)
550 nonRegionCells.append(cellI);
553 return nonRegionCells.shrink();
559 void getInterfaceSizes
563 const bool sumParallel,
575 label ownRegion = cellRegion[mesh.
faceOwner()[faceI]];
578 if (ownRegion != neiRegion)
582 min(ownRegion, neiRegion),
583 max(ownRegion, neiRegion)
588 if (iter != interfaceSizes.
end())
608 coupledRegion[i] = cellRegion[cellI];
615 label ownRegion = cellRegion[mesh.
faceOwner()[faceI]];
616 label neiRegion = coupledRegion[i];
618 if (ownRegion != neiRegion)
622 min(ownRegion, neiRegion),
623 max(ownRegion, neiRegion)
628 if (iter != interfaceSizes.
end())
659 interfaceSizes.
find(slaveIter.key());
661 if (masterIter != interfaceSizes.
end())
663 masterIter() += slaveIter();
667 interfaceSizes.
insert(slaveIter.key(), slaveIter());
682 toSlave << interfaceSizes;
690 toMaster << interfaceSizes;
695 fromMaster >> interfaceSizes;
701 interfaces = interfaceSizes.
toc();
733 Info<<
"Testing:" << io.objectPath() <<
endl;
738 Info<<
"Writing dummy " << regionName/io.name() <<
endl;
741 dummyDict.
add(
"divSchemes", divDict);
743 dummyDict.
add(
"gradSchemes", gradDict);
745 dummyDict.
add(
"laplacianSchemes", laplDict);
754 mesh.time().system(),
765 Info<<
"Writing dummy " << regionName/io.name() <<
endl;
773 labelList coupledRegion(mesh.nFaces()-mesh.nInternalFaces());
777 label cellI = mesh.faceOwner()[i+mesh.nInternalFaces()];
778 coupledRegion[i] = cellRegion[cellI];
790 labelList cellsToRemove(getNonRegionCells(cellRegion, regionI));
795 labelList exposedFaces = cellRemover.getExposedFaces(cellsToRemove);
800 label faceI = exposedFaces[i];
802 label ownRegion = cellRegion[mesh.faceOwner()[faceI]];
803 label neiRegion = -1;
805 if (mesh.isInternalFace(faceI))
807 neiRegion = cellRegion[mesh.faceNeighbour()[faceI]];
811 neiRegion = coupledRegion[faceI-mesh.nInternalFaces()];
814 label otherRegion = -1;
816 if (ownRegion == regionI && neiRegion != regionI)
818 otherRegion = neiRegion;
820 else if (ownRegion != regionI && neiRegion == regionI)
822 otherRegion = ownRegion;
827 <<
"Exposed face:" << faceI
828 <<
" fc:" << mesh.faceCentres()[faceI]
829 <<
" has owner region " << ownRegion
830 <<
" and neighbour region " << neiRegion
831 <<
" when handling region:" << regionI
835 if (otherRegion != -1)
840 if (regionI < otherRegion)
842 exposedPatchIDs[i] = interfaceToPatch[
interface];
846 exposedPatchIDs[i] = interfaceToPatch[
interface]+1;
852 cellRemover.setRefinement
866 mesh.time().timeName(),
878 void createAndWriteRegion
885 const word& newMeshInstance
888 Info<<
"Creating mesh for region " << regionI
889 <<
' ' << regionNames[regionI] <<
endl;
898 regionNames[regionI],
907 addedPatches.
insert(iter());
908 addedPatches.
insert(iter()+1);
914 newMesh().updateMesh(map());
917 subsetVolFields<volScalarField>
925 subsetVolFields<volVectorField>
933 subsetVolFields<volSphericalTensorField>
941 subsetVolFields<volSymmTensorField>
949 subsetVolFields<volTensorField>
958 subsetSurfaceFields<surfaceScalarField>
965 subsetSurfaceFields<surfaceVectorField>
972 subsetSurfaceFields<surfaceSphericalTensorField>
979 subsetSurfaceFields<surfaceSymmTensorField>
986 subsetSurfaceFields<surfaceTensorField>
996 newPatches.checkParallelSync(
true);
1002 labelList oldToNew(newPatches.size(), -1);
1005 Info<<
"Deleting empty patches" <<
endl;
1008 forAll(newPatches, patchI)
1010 const polyPatch& pp = newPatches[patchI];
1012 if (!isA<processorPolyPatch>(pp))
1016 oldToNew[patchI] = newI++;
1022 forAll(newPatches, patchI)
1024 const polyPatch& pp = newPatches[patchI];
1026 if (isA<processorPolyPatch>(pp) && pp.size())
1028 oldToNew[patchI] = newI++;
1032 const label nNewPatches = newI;
1037 if (oldToNew[patchI] == -1)
1039 oldToNew[patchI] = newI++;
1043 reorderPatches(newMesh(), oldToNew, nNewPatches);
1048 newMesh().setInstance(newMeshInstance);
1052 Info<<
"Writing addressing to base mesh" <<
endl;
1058 "pointRegionAddressing",
1059 newMesh().facesInstance(),
1060 newMesh().meshSubDir,
1068 Info<<
"Writing map " << pointProcAddressing.name()
1069 <<
" from region" << regionI
1070 <<
" points back to base mesh." <<
endl;
1071 pointProcAddressing.
write();
1077 "faceRegionAddressing",
1078 newMesh().facesInstance(),
1079 newMesh().meshSubDir,
1087 forAll(faceProcAddressing, faceI)
1091 label oldFaceI = map().faceMap()[faceI];
1095 map().cellMap()[newMesh().faceOwner()[faceI]]
1096 == mesh.faceOwner()[oldFaceI]
1099 faceProcAddressing[faceI] = oldFaceI+1;
1103 faceProcAddressing[faceI] = -(oldFaceI+1);
1106 Info<<
"Writing map " << faceProcAddressing.name()
1107 <<
" from region" << regionI
1108 <<
" faces back to base mesh." <<
endl;
1109 faceProcAddressing.
write();
1115 "cellRegionAddressing",
1116 newMesh().facesInstance(),
1117 newMesh().meshSubDir,
1125 Info<<
"Writing map " <<cellProcAddressing.name()
1126 <<
" from region" << regionI
1127 <<
" cells back to base mesh." <<
endl;
1128 cellProcAddressing.
write();
1140 const label nCellRegions,
1153 forAll(interfaces, interI)
1155 const edge&
e = interfaces[interI];
1157 if (interfaceSizes[e] > 0)
1159 const word inter1 = regionNames[e[0]] +
"_to_" + regionNames[e[1]];
1160 const word inter2 = regionNames[e[1]] +
"_to_" + regionNames[e[0]];
1175 label patchI = addPatch(mesh, patch1);
1189 addPatch(mesh, patch2);
1191 Info<<
"For interface between region " << e[0]
1192 <<
" and " << e[1] <<
" added patch " << patchI
1196 interfaceToPatch.
insert(e, patchI);
1199 return interfaceToPatch;
1204 label findCorrespondingRegion
1208 const label nCellRegions,
1210 const label minOverlapSize
1216 forAll(cellRegion, cellI)
1218 if (existingZoneID[cellI] == zoneI)
1220 cellsInZone[cellRegion[cellI]]++;
1228 label regionI =
findMax(cellsInZone);
1231 if (cellsInZone[regionI] < minOverlapSize)
1239 forAll(cellRegion, cellI)
1241 if (cellRegion[cellI] == regionI && existingZoneID[cellI] != zoneI)
1340 const cellZone& cz = cellZones[zoneI];
1344 label cellI = cz[i];
1345 if (zoneID[cellI] == -1)
1347 zoneID[cellI] = zoneI;
1352 <<
"Cell " << cellI <<
" with cell centre "
1354 <<
" is multiple zones. This is not allowed." << endl
1355 <<
"It is in zone " << cellZones[zoneID[cellI]].
name()
1356 <<
" and in zone " << cellZones[zoneI].
name()
1375 const bool sloppyCellZones,
1378 const label nCellRegions,
1388 regionToZone.
setSize(nCellRegions, -1);
1389 regionNames.
setSize(nCellRegions);
1395 getZoneID(mesh, cellZones, zoneID, neiZoneID);
1407 if (zoneNames[procI] != zoneNames[0])
1410 <<
"cellZones not synchronised across processors." << endl
1411 <<
"Master has cellZones " << zoneNames[0] << endl
1412 <<
"Processor " << procI
1413 <<
" has cellZones " << zoneNames[procI]
1422 cellZones[zoneI].size(),
1429 if (sloppyCellZones)
1431 Info<<
"Trying to match regions to existing cell zones;"
1432 <<
" region can be subset of cell zone." <<
nl <<
endl;
1436 label regionI = findCorrespondingRegion
1442 label(0.5*zoneSizes[zoneI])
1447 Info<<
"Sloppily matched region " << regionI
1449 <<
" to zone " << zoneI <<
" size " << zoneSizes[zoneI]
1451 zoneToRegion[zoneI] = regionI;
1452 regionToZone[regionI] = zoneI;
1453 regionNames[regionI] = cellZones[zoneI].
name();
1459 Info<<
"Trying to match regions to existing cell zones." <<
nl <<
endl;
1463 label regionI = findCorrespondingRegion
1474 zoneToRegion[zoneI] = regionI;
1475 regionToZone[regionI] = zoneI;
1476 regionNames[regionI] = cellZones[zoneI].
name();
1481 forAll(regionToZone, regionI)
1483 if (regionToZone[regionI] == -1)
1485 regionNames[regionI] =
"domain" +
Foam::name(regionI);
1491 void writeCellToRegion(
const fvMesh& mesh,
const labelList& cellRegion)
1508 cellToRegion.write();
1510 Info<<
"Writing region per cell file (for manual decomposition) to "
1511 << cellToRegion.objectPath() <<
nl <<
endl;
1528 zeroGradientFvPatchScalarField::typeName
1530 forAll(cellRegion, cellI)
1532 cellToRegion[cellI] = cellRegion[cellI];
1534 cellToRegion.write();
1536 Info<<
"Writing region per cell as volScalarField to "
1537 << cellToRegion.objectPath() <<
nl <<
endl;
1545 int main(
int argc,
char *argv[])
1560 runTime.functionObjects().off();
1564 word blockedFacesName;
1567 blockedFacesName =
args.
option(
"blockedFaces");
1568 Info<<
"Reading blocked internal faces from faceSet "
1569 << blockedFacesName <<
nl <<
endl;
1584 (useCellZonesOnly || useCellZonesFile)
1592 <<
"You cannot specify both -cellZonesOnly or -cellZonesFileOnly"
1593 <<
" (which specify complete split)"
1594 <<
" in combination with -blockedFaces or -cellZones"
1595 <<
" (which imply a split based on topology)"
1601 if (insidePoint && largestOnly)
1604 <<
"You cannot specify both -largestOnly"
1605 <<
" (keep region with most cells)"
1606 <<
" and -insidePoint (keep region containing point)"
1617 getZoneID(mesh, cellZones, zoneID, neiZoneID);
1633 label nCellRegions = 0;
1634 if (useCellZonesOnly)
1636 Info<<
"Using current cellZones to split mesh into regions."
1637 <<
" This requires all"
1638 <<
" cells to be in one and only one cellZone." <<
nl <<
endl;
1640 label unzonedCellI =
findIndex(zoneID, -1);
1641 if (unzonedCellI != -1)
1644 <<
"For the cellZonesOnly option all cells "
1645 <<
"have to be in a cellZone." << endl
1646 <<
"Cell " << unzonedCellI
1648 <<
" is not in a cellZone. There might be more unzoned cells."
1651 cellRegion = zoneID;
1652 nCellRegions =
gMax(cellRegion)+1;
1653 regionToZone.
setSize(nCellRegions);
1654 regionNames.
setSize(nCellRegions);
1656 for (label regionI = 0; regionI < nCellRegions; regionI++)
1658 regionToZone[regionI] = regionI;
1659 zoneToRegion[regionI] = regionI;
1660 regionNames[regionI] = cellZones[regionI].
name();
1663 else if (useCellZonesFile)
1666 Info<<
"Reading split from cellZones file " << zoneFile << endl
1667 <<
"This requires all"
1668 <<
" cells to be in one and only one cellZone." <<
nl <<
endl;
1687 getZoneID(mesh, newCellZones, newZoneID, newNeiZoneID);
1689 label unzonedCellI =
findIndex(newZoneID, -1);
1690 if (unzonedCellI != -1)
1693 <<
"For the cellZonesFileOnly option all cells "
1694 <<
"have to be in a cellZone." << endl
1695 <<
"Cell " << unzonedCellI
1697 <<
" is not in a cellZone. There might be more unzoned cells."
1700 cellRegion = newZoneID;
1701 nCellRegions =
gMax(cellRegion)+1;
1702 zoneToRegion.
setSize(newCellZones.size(), -1);
1703 regionToZone.
setSize(nCellRegions);
1704 regionNames.
setSize(nCellRegions);
1705 for (label regionI = 0; regionI < nCellRegions; regionI++)
1707 regionToZone[regionI] = regionI;
1708 zoneToRegion[regionI] = regionI;
1709 regionNames[regionI] = newCellZones[regionI].name();
1721 if (blockedFacesName.size())
1723 faceSet blockedFaceSet(mesh, blockedFacesName);
1726 <<
" blocked faces from set " << blockedFacesName <<
nl <<
endl;
1732 blockedFace[iter.key()] =
true;
1746 if (zoneID[own] != zoneID[nei])
1748 blockedFace[faceI] =
true;
1757 if (zoneID[mesh.
faceOwner()[faceI]] != neiZoneID[i])
1759 blockedFace[faceI] =
true;
1766 nCellRegions = regions.nRegions();
1783 Info<< endl <<
"Number of regions:" << nCellRegions <<
nl <<
endl;
1787 writeCellToRegion(mesh, cellRegion);
1796 forAll(cellRegion, cellI)
1798 regionSizes[cellRegion[cellI]]++;
1800 forAll(regionSizes, regionI)
1805 Info<<
"Region\tCells" <<
nl
1806 <<
"------\t-----" <<
endl;
1808 forAll(regionSizes, regionI)
1810 Info<< regionI <<
'\t' << regionSizes[regionI] <<
nl;
1817 Info<<
"Region\tZone\tName" <<
nl
1818 <<
"------\t----\t----" <<
endl;
1819 forAll(regionToZone, regionI)
1821 Info<< regionI <<
'\t' << regionToZone[regionI] <<
'\t'
1822 << regionNames[regionI] <<
nl;
1847 Info<<
"Sizes inbetween regions:" <<
nl <<
nl
1848 <<
"Region\tRegion\tFaces" <<
nl
1849 <<
"------\t------\t-----" <<
endl;
1851 forAll(interfaces, interI)
1853 const edge& e = interfaces[interI];
1855 Info<< e[0] <<
'\t' << e[1] <<
'\t' << interfaceSizes[
e] <<
nl;
1912 if (nCellRegions == 1)
1914 Info<<
"Only one region. Doing nothing." <<
endl;
1916 else if (makeCellZones)
1918 Info<<
"Putting cells into cellZones instead of splitting mesh."
1923 for (label regionI = 0; regionI < nCellRegions; regionI++)
1925 label zoneI = regionToZone[regionI];
1929 Info<<
" Region " << regionI <<
" : corresponds to existing"
1931 << zoneI <<
' ' << cellZones[zoneI].
name() <<
endl;
1944 zoneI = cellZones.
size();
1963 Info<<
" Region " << regionI <<
" : created new cellZone "
1964 << zoneI <<
' ' << cellZones[zoneI].
name() <<
endl;
1979 Info<<
"Writing cellZones as new mesh to time " << runTime.timeName()
1988 Info<<
"Writing cellSets corresponding to cellZones." <<
nl <<
endl;
1992 const cellZone& cz = cellZones[zoneI];
2034 label cellI = mesh.
findCell(insidePoint);
2036 Info<<
nl <<
"Found point " << insidePoint <<
" in cell " << cellI
2041 regionI = cellRegion[cellI];
2047 <<
"Subsetting region " << regionI
2048 <<
" containing point " << insidePoint <<
endl;
2053 <<
"Point " << insidePoint
2054 <<
" is not inside the mesh." <<
nl
2055 <<
"Bounding box of the mesh:" << mesh.
bounds()
2059 createAndWriteRegion
2066 (overwrite ? oldInstance : runTime.timeName())
2069 else if (largestOnly)
2071 label regionI =
findMax(regionSizes);
2074 <<
"Subsetting region " << regionI
2075 <<
" of size " << regionSizes[regionI] <<
endl;
2077 createAndWriteRegion
2084 (overwrite ? oldInstance : runTime.timeName())
2090 for (label regionI = 0; regionI < nCellRegions; regionI++)
2093 <<
"Region " << regionI <<
nl
2094 <<
"-------- " <<
endl;
2096 createAndWriteRegion
2103 (overwrite ? oldInstance : runTime.timeName())