zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: Added option for explicitly providing regions. Restructured code.

Browse Source
This commit is contained in:
mattijs
2010-04-01 14:20:27 +01:00
parent 39ea47af82
commit fd7e9c9ae4
1 changed files with 376 additions and 292 deletions
Download Patch File Download Diff File Expand all files Collapse all files

668
applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
View File

@ -39,17 +39,29 @@ Description
- mesh with cells put into cellZones (-makeCellZones)
Note:
- cellZonesOnly does not do a walk and uses the cellZones only. Use
this if you don't mind having disconnected domains in a single region.
This option requires all cells to be in one (and one only) cellZone.
- cellZonesFileOnly behaves like -cellZonesOnly but reads the cellZones
from the specified file. This allows one to explicitly specify the region
distribution and still have multiple cellZones per region.
- useCellZonesOnly does not do a walk and uses the cellZones only. Use
this if you don't mind having disconnected domains in a single region.
This option requires all cells to be in one (and one only) cellZone.
- Should work in parallel.
cellZones can differ on either side of processor boundaries in which case
the faces get moved from processor patch to directMapped patch. Not
ery well tested.
very well tested.
- If a cell zone gets split into more than one region it can detect
the largest matching region (-sloppyCellZones). This will accept any
region that covers more than 50% of the zone. It has to be a subset
so cannot have any cells in any other zone.
- useCellZonesOnly does not do a walk and uses the cellZones only. Use
this if you don't mind having disconnected domains in a single region.
This option requires all cells to be in one (and one only) cellZone.
- writes maps like decomposePar back to original mesh:
- pointRegionAddressing : for every point in this region the point in
the original mesh
@ -1137,63 +1149,6 @@ EdgeMap<label> addRegionPatches
}
//// Checks if regionI in cellRegion is subset of existing cellZone. Returns -1
//// if no zone found, zone otherwise
//label findCorrespondingSubZone
//(
// const cellZoneMesh& cellZones,
// const labelList& existingZoneID,
// const labelList& cellRegion,
// const label regionI
//)
//{
// // Zone corresponding to region. No corresponding zone.
// label zoneI = labelMax;
//
// labelList regionCells = findIndices(cellRegion, regionI);
//
// if (regionCells.empty())
// {
// // My local portion is empty. Maps to any empty cellZone. Mark with
// // special value which can get overwritten by other processors.
// zoneI = -1;
// }
// else
// {
// // Get zone for first element.
// zoneI = existingZoneID[regionCells[0]];
//
// if (zoneI == -1)
// {
// zoneI = labelMax;
// }
// else
// {
// // 1. All regionCells in zoneI?
// forAll(regionCells, i)
// {
// if (existingZoneID[regionCells[i]] != zoneI)
// {
// zoneI = labelMax;
// break;
// }
// }
// }
// }
//
// // Determine same zone over all processors.
// reduce(zoneI, maxOp<label>());
//
// if (zoneI == labelMax)
// {
// // Cells in region that are not in zoneI
// zoneI = -1;
// }
//
// return zoneI;
//}
// Find region that covers most of cell zone
label findCorrespondingRegion
(
@ -1314,62 +1269,20 @@ label findCorrespondingRegion
//}
// Main program:
int main(int argc, char *argv[])
// Get zone per cell
// - non-unique zoning
// - coupled zones
void getZoneID
(
const polyMesh& mesh,
const cellZoneMesh& cellZones,
labelList& zoneID,
labelList& neiZoneID
)
{
# include "addOverwriteOption.H"
argList::addBoolOption("cellZones");
argList::addBoolOption("cellZonesOnly");
argList::addOption("blockedFaces", "faceSet");
argList::addBoolOption("makeCellZones");
argList::addBoolOption("largestOnly");
argList::addOption("insidePoint", "point");
argList::addBoolOption("detectOnly");
argList::addBoolOption("sloppyCellZones");
# include "setRootCase.H"
# include "createTime.H"
runTime.functionObjects().off();
# include "createMesh.H"
const word oldInstance = mesh.pointsInstance();
word blockedFacesName;
if (args.optionReadIfPresent("blockedFaces", blockedFacesName))
{
Info<< "Reading blocked internal faces from faceSet "
<< blockedFacesName << nl << endl;
}
const bool makeCellZones = args.optionFound("makeCellZones");
const bool largestOnly = args.optionFound("largestOnly");
const bool insidePoint = args.optionFound("insidePoint");
const bool useCellZones = args.optionFound("cellZones");
const bool useCellZonesOnly = args.optionFound("cellZonesOnly");
const bool overwrite = args.optionFound("overwrite");
const bool detectOnly = args.optionFound("detectOnly");
const bool sloppyCellZones = args.optionFound("sloppyCellZones");
if (insidePoint && largestOnly)
{
FatalErrorIn(args.executable())
<< "You cannot specify both -largestOnly"
<< " (keep region with most cells)"
<< " and -insidePoint (keep region containing point)"
<< exit(FatalError);
}
const cellZoneMesh& cellZones = mesh.cellZones();
// Collect zone per cell
// ~~~~~~~~~~~~~~~~~~~~~
// - non-unique zoning
// - coupled zones
// Existing zoneID
labelList zoneID(mesh.nCells(), -1);
zoneID.setSize(mesh.nCells());
zoneID = -1;
forAll(cellZones, zoneI)
{
@ -1384,7 +1297,7 @@ int main(int argc, char *argv[])
}
else
{
FatalErrorIn(args.executable())
FatalErrorIn("getZoneID(..)")
<< "Cell " << cellI << " with cell centre "
<< mesh.cellCentres()[cellI]
<< " is multiple zones. This is not allowed." << endl
@ -1396,175 +1309,42 @@ int main(int argc, char *argv[])
}
// Neighbour zoneID.
labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
neiZoneID.setSize(mesh.nFaces()-mesh.nInternalFaces());
forAll(neiZoneID, i)
{
neiZoneID[i] = zoneID[mesh.faceOwner()[i+mesh.nInternalFaces()]];
}
syncTools::swapBoundaryFaceList(mesh, neiZoneID, false);
}
// Determine connected regions
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
void matchRegions
(
const bool sloppyCellZones,
const polyMesh& mesh,
// Mark additional faces that are blocked
boolList blockedFace;
const label nCellRegions,
const labelList& cellRegion,
// Read from faceSet
if (blockedFacesName.size())
{
faceSet blockedFaceSet(mesh, blockedFacesName);
Info<< "Read " << returnReduce(blockedFaceSet.size(), sumOp<label>())
<< " blocked faces from set " << blockedFacesName << nl << endl;
labelList& regionToZone,
wordList& regionNames,
labelList& zoneToRegion
)
{
const cellZoneMesh& cellZones = mesh.cellZones();
blockedFace.setSize(mesh.nFaces(), false);
forAllConstIter(faceSet, blockedFaceSet, iter)
{
blockedFace[iter.key()] = true;
}
}
// Imply from differing cellZones
if (useCellZones)
{
blockedFace.setSize(mesh.nFaces(), false);
for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
{
label own = mesh.faceOwner()[faceI];
label nei = mesh.faceNeighbour()[faceI];
if (zoneID[own] != zoneID[nei])
{
blockedFace[faceI] = true;
}
}
// Different cellZones on either side of processor patch.
forAll(neiZoneID, i)
{
label faceI = i+mesh.nInternalFaces();
if (zoneID[mesh.faceOwner()[faceI]] != neiZoneID[i])
{
blockedFace[faceI] = true;
}
}
}
// Determine per cell the region it belongs to
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// cellRegion is the labelList with the region per cell.
labelList cellRegion;
label nCellRegions = 0;
if (useCellZonesOnly)
{
label unzonedCellI = findIndex(zoneID, -1);
if (unzonedCellI != -1)
{
FatalErrorIn(args.executable())
<< "For the cellZonesOnly option all cells "
<< "have to be in a cellZone." << endl
<< "Cell " << unzonedCellI
<< " at" << mesh.cellCentres()[unzonedCellI]
<< " is not in a cellZone. There might be more unzoned cells."
<< exit(FatalError);
}
cellRegion = zoneID;
nCellRegions = gMax(cellRegion)+1;
}
else
{
// Do a topological walk to determine regions
regionSplit regions(mesh, blockedFace);
nCellRegions = regions.nRegions();
cellRegion.transfer(regions);
}
Info<< endl << "Number of regions:" << nCellRegions << nl << endl;
// Write to manual decomposition option
{
labelIOList cellToRegion
(
IOobject
(
"cellToRegion",
mesh.facesInstance(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
cellRegion
);
cellToRegion.write();
Info<< "Writing region per cell file (for manual decomposition) to "
<< cellToRegion.objectPath() << nl << endl;
}
// Write for postprocessing
{
volScalarField cellToRegion
(
IOobject
(
"cellToRegion",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh,
dimensionedScalar("zero", dimless, 0),
zeroGradientFvPatchScalarField::typeName
);
forAll(cellRegion, cellI)
{
cellToRegion[cellI] = cellRegion[cellI];
}
cellToRegion.write();
Info<< "Writing region per cell as volScalarField to "
<< cellToRegion.objectPath() << nl << endl;
}
// Sizes per region
// ~~~~~~~~~~~~~~~~
labelList regionSizes(nCellRegions, 0);
forAll(cellRegion, cellI)
{
regionSizes[cellRegion[cellI]]++;
}
forAll(regionSizes, regionI)
{
reduce(regionSizes[regionI], sumOp<label>());
}
Info<< "Region\tCells" << nl
<< "------\t-----" << endl;
forAll(regionSizes, regionI)
{
Info<< regionI << '\t' << regionSizes[regionI] << nl;
}
Info<< endl;
regionToZone.setSize(nCellRegions, -1);
regionNames.setSize(nCellRegions);
zoneToRegion.setSize(cellZones.size(), -1);
// Get current per cell zoneID
labelList zoneID(mesh.nCells(), -1);
labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
getZoneID(mesh, cellZones, zoneID, neiZoneID);
// Sizes per cellzone
// ~~~~~~~~~~~~~~~~~~
labelList zoneSizes(cellZones.size(), 0);
if (useCellZones || makeCellZones || sloppyCellZones)
{
List<wordList> zoneNames(Pstream::nProcs());
zoneNames[Pstream::myProcNo()] = cellZones.names();
@ -1575,7 +1355,7 @@ int main(int argc, char *argv[])
{
if (zoneNames[procI] != zoneNames[0])
{
FatalErrorIn(args.executable())
FatalErrorIn("matchRegions(..)")
<< "cellZones not synchronised across processors." << endl
<< "Master has cellZones " << zoneNames[0] << endl
<< "Processor " << procI
@ -1595,16 +1375,6 @@ int main(int argc, char *argv[])
}
// Whether region corresponds to a cellzone
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Region per zone
labelList regionToZone(nCellRegions, -1);
// Name of region
wordList regionNames(nCellRegions);
// Zone to region
labelList zoneToRegion(cellZones.size(), -1);
if (sloppyCellZones)
{
Info<< "Trying to match regions to existing cell zones;"
@ -1624,7 +1394,7 @@ int main(int argc, char *argv[])
if (regionI != -1)
{
Info<< "Sloppily matched region " << regionI
<< " size " << regionSizes[regionI]
//<< " size " << regionSizes[regionI]
<< " to zone " << zoneI << " size " << zoneSizes[zoneI]
<< endl;
zoneToRegion[zoneI] = regionI;
@ -1664,6 +1434,330 @@ int main(int argc, char *argv[])
regionNames[regionI] = "domain" + Foam::name(regionI);
}
}
}
void writeCellToRegion(const fvMesh& mesh, const labelList& cellRegion)
{
// Write to manual decomposition option
{
labelIOList cellToRegion
(
IOobject
(
"cellToRegion",
mesh.facesInstance(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
cellRegion
);
cellToRegion.write();
Info<< "Writing region per cell file (for manual decomposition) to "
<< cellToRegion.objectPath() << nl << endl;
}
// Write for postprocessing
{
volScalarField cellToRegion
(
IOobject
(
"cellToRegion",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh,
dimensionedScalar("zero", dimless, 0),
zeroGradientFvPatchScalarField::typeName
);
forAll(cellRegion, cellI)
{
cellToRegion[cellI] = cellRegion[cellI];
}
cellToRegion.write();
Info<< "Writing region per cell as volScalarField to "
<< cellToRegion.objectPath() << nl << endl;
}
}
// Main program:
int main(int argc, char *argv[])
{
# include "addOverwriteOption.H"
argList::addBoolOption("cellZones");
argList::addBoolOption("cellZonesOnly");
argList::addOption("cellZonesFileOnly", "cellZonesName");
argList::addOption("blockedFaces", "faceSet");
argList::addBoolOption("makeCellZones");
argList::addBoolOption("largestOnly");
argList::addOption("insidePoint", "point");
argList::addBoolOption("detectOnly");
argList::addBoolOption("sloppyCellZones");
# include "setRootCase.H"
# include "createTime.H"
runTime.functionObjects().off();
# include "createMesh.H"
const word oldInstance = mesh.pointsInstance();
word blockedFacesName;
if (args.optionReadIfPresent("blockedFaces", blockedFacesName))
{
Info<< "Reading blocked internal faces from faceSet "
<< blockedFacesName << nl << endl;
}
const bool makeCellZones = args.optionFound("makeCellZones");
const bool largestOnly = args.optionFound("largestOnly");
const bool insidePoint = args.optionFound("insidePoint");
const bool useCellZones = args.optionFound("cellZones");
const bool useCellZonesOnly = args.optionFound("cellZonesOnly");
const bool useCellZonesFile = args.optionFound("cellZonesFileOnly");
const bool overwrite = args.optionFound("overwrite");
const bool detectOnly = args.optionFound("detectOnly");
const bool sloppyCellZones = args.optionFound("sloppyCellZones");
if
(
(useCellZonesOnly || useCellZonesFile)
&& (
blockedFacesName != word::null
|| useCellZones
)
)
{
FatalErrorIn(args.executable())
<< "You cannot specify both -cellZonesOnly or -cellZonesFileOnly"
<< " (which specify complete split)"
<< " in combination with -blockedFaces or -cellZones"
<< " (which imply a split based on topology)"
<< exit(FatalError);
}
if (insidePoint && largestOnly)
{
FatalErrorIn(args.executable())
<< "You cannot specify both -largestOnly"
<< " (keep region with most cells)"
<< " and -insidePoint (keep region containing point)"
<< exit(FatalError);
}
const cellZoneMesh& cellZones = mesh.cellZones();
// Existing zoneID
labelList zoneID(mesh.nCells(), -1);
// Neighbour zoneID.
labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
getZoneID(mesh, cellZones, zoneID, neiZoneID);
// Determine per cell the region it belongs to
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// cellRegion is the labelList with the region per cell.
labelList cellRegion;
// Region per zone
labelList regionToZone;
// Name of region
wordList regionNames;
// Zone to region
labelList zoneToRegion;
label nCellRegions = 0;
if (useCellZonesOnly)
{
Info<< "Using current cellZones to split mesh into regions."
<< " This requires all"
<< " cells to be in one and only one cellZone." << nl << endl;
label unzonedCellI = findIndex(zoneID, -1);
if (unzonedCellI != -1)
{
FatalErrorIn(args.executable())
<< "For the cellZonesOnly option all cells "
<< "have to be in a cellZone." << endl
<< "Cell " << unzonedCellI
<< " at" << mesh.cellCentres()[unzonedCellI]
<< " is not in a cellZone. There might be more unzoned cells."
<< exit(FatalError);
}
cellRegion = zoneID;
nCellRegions = gMax(cellRegion)+1;
regionToZone.setSize(nCellRegions);
regionNames.setSize(nCellRegions);
zoneToRegion.setSize(cellZones.size(), -1);
for (label regionI = 0; regionI < nCellRegions; regionI++)
{
regionToZone[regionI] = regionI;
zoneToRegion[regionI] = regionI;
regionNames[regionI] = cellZones[regionI].name();
}
}
else if (useCellZonesFile)
{
const word zoneFile = args.option("cellZonesFileOnly");
Info<< "Reading split from cellZones file " << zoneFile << endl
<< "This requires all"
<< " cells to be in one and only one cellZone." << nl << endl;
cellZoneMesh newCellZones
(
IOobject
(
zoneFile,
mesh.facesInstance(),
polyMesh::meshSubDir,
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
),
mesh
);
labelList newZoneID(mesh.nCells(), -1);
labelList newNeiZoneID(mesh.nFaces()-mesh.nInternalFaces());
getZoneID(mesh, newCellZones, newZoneID, newNeiZoneID);
label unzonedCellI = findIndex(newZoneID, -1);
if (unzonedCellI != -1)
{
FatalErrorIn(args.executable())
<< "For the cellZonesFileOnly option all cells "
<< "have to be in a cellZone." << endl
<< "Cell " << unzonedCellI
<< " at" << mesh.cellCentres()[unzonedCellI]
<< " is not in a cellZone. There might be more unzoned cells."
<< exit(FatalError);
}
cellRegion = newZoneID;
nCellRegions = gMax(cellRegion)+1;
zoneToRegion.setSize(newCellZones.size(), -1);
regionToZone.setSize(nCellRegions);
regionNames.setSize(nCellRegions);
for (label regionI = 0; regionI < nCellRegions; regionI++)
{
regionToZone[regionI] = regionI;
zoneToRegion[regionI] = regionI;
regionNames[regionI] = newCellZones[regionI].name();
}
}
else
{
// Determine connected regions
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Mark additional faces that are blocked
boolList blockedFace;
// Read from faceSet
if (blockedFacesName.size())
{
faceSet blockedFaceSet(mesh, blockedFacesName);
Info<< "Read "
<< returnReduce(blockedFaceSet.size(), sumOp<label>())
<< " blocked faces from set " << blockedFacesName << nl << endl;
blockedFace.setSize(mesh.nFaces(), false);
forAllConstIter(faceSet, blockedFaceSet, iter)
{
blockedFace[iter.key()] = true;
}
}
// Imply from differing cellZones
if (useCellZones)
{
blockedFace.setSize(mesh.nFaces(), false);
for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
{
label own = mesh.faceOwner()[faceI];
label nei = mesh.faceNeighbour()[faceI];
if (zoneID[own] != zoneID[nei])
{
blockedFace[faceI] = true;
}
}
// Different cellZones on either side of processor patch.
forAll(neiZoneID, i)
{
label faceI = i+mesh.nInternalFaces();
if (zoneID[mesh.faceOwner()[faceI]] != neiZoneID[i])
{
blockedFace[faceI] = true;
}
}
}
// Do a topological walk to determine regions
regionSplit regions(mesh, blockedFace);
nCellRegions = regions.nRegions();
cellRegion.transfer(regions);
// Make up region names. If possible match them to existing zones.
matchRegions
(
sloppyCellZones,
mesh,
nCellRegions,
cellRegion,
regionToZone,
regionNames,
zoneToRegion
);
}
Info<< endl << "Number of regions:" << nCellRegions << nl << endl;
// Write decomposition to file
writeCellToRegion(mesh, cellRegion);
// Sizes per region
// ~~~~~~~~~~~~~~~~
labelList regionSizes(nCellRegions, 0);
forAll(cellRegion, cellI)
{
regionSizes[cellRegion[cellI]]++;
}
forAll(regionSizes, regionI)
{
reduce(regionSizes[regionI], sumOp<label>());
}
Info<< "Region\tCells" << nl
<< "------\t-----" << endl;
forAll(regionSizes, regionI)
{
Info<< regionI << '\t' << regionSizes[regionI] << nl;
}
Info<< endl;
// Print region to zone
@ -1676,16 +1770,6 @@ int main(int argc, char *argv[])
}
Info<< endl;
//// Print zone to region
//Info<< "Zone\tName\tRegion" << nl
// << "----\t----\t------" << endl;
//forAll(zoneToRegion, zoneI)
//{
// Info<< zoneI << '\t' << cellZones[zoneI].name() << '\t'
// << zoneToRegion[zoneI] << nl;
//}
//Info<< endl;
// Since we're going to mess with patches make sure all non-processor ones