2205 lines
59 KiB
C
2205 lines
59 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | www.openfoam.com
|
|
\\/ M anipulation |
|
|
-------------------------------------------------------------------------------
|
|
Copyright (C) 2011-2017 OpenFOAM Foundation
|
|
Copyright (C) 2015-2021 OpenCFD Ltd.
|
|
-------------------------------------------------------------------------------
|
|
License
|
|
This file is part of OpenFOAM.
|
|
|
|
OpenFOAM is free software: you can redistribute it and/or modify it
|
|
under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
|
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
Application
|
|
splitMeshRegions
|
|
|
|
Group
|
|
grpMeshManipulationUtilities
|
|
|
|
Description
|
|
Splits mesh into multiple regions.
|
|
|
|
Each region is defined as a domain whose cells can all be reached by
|
|
cell-face-cell walking without crossing
|
|
- boundary faces
|
|
- additional faces from faceset (-blockedFaces faceSet).
|
|
- any face between differing cellZones (-cellZones)
|
|
|
|
Output is:
|
|
- volScalarField with regions as different scalars (-detectOnly)
|
|
or
|
|
- mesh with multiple regions and mapped patches. These patches
|
|
either cover the whole interface between two region (default) or
|
|
only part according to faceZones (-useFaceZones)
|
|
or
|
|
- mesh with cells put into cellZones (-makeCellZones)
|
|
|
|
Note:
|
|
|
|
- multiple cellZones can be combined into a single cellZone (cluster)
|
|
for further analysis using the 'combineZones' option:
|
|
-combineZones '((zoneA "zoneB.*")(none otherZone))
|
|
This can be combined with e.g. 'cellZones' or 'cellZonesOnly'. The
|
|
corresponding region name will be the names of the cellZones combined e.g.
|
|
zoneA_zoneB0_zoneB1
|
|
|
|
- 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.
|
|
|
|
- prefixRegion prefixes all normal patches with region name (interface
|
|
(patches already have region name prefix)
|
|
|
|
- Should work in parallel.
|
|
cellZones can differ on either side of processor boundaries in which case
|
|
the faces get moved from processor patch to mapped patch. Not 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.
|
|
|
|
- If explicitly a single region has been selected (-largestOnly or
|
|
-insidePoint) its region name will be either
|
|
- name of a cellZone it matches to or
|
|
- "largestOnly" respectively "insidePoint" or
|
|
- polyMesh::defaultRegion if additionally -overwrite
|
|
(so it will overwrite the input mesh!)
|
|
|
|
- writes maps like decomposePar back to original mesh:
|
|
- pointRegionAddressing : for every point in this region the point in
|
|
the original mesh
|
|
- cellRegionAddressing : ,, cell ,, cell ,,
|
|
- faceRegionAddressing : ,, face ,, face in
|
|
the original mesh + 'turning index'. For a face in the same orientation
|
|
this is the original facelabel+1, for a turned face this is -facelabel-1
|
|
- boundaryRegionAddressing : for every patch in this region the
|
|
patch in the original mesh (or -1 if added patch)
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "SortableList.H"
|
|
#include "argList.H"
|
|
#include "regionSplit.H"
|
|
#include "fvMeshSubset.H"
|
|
#include "IOobjectList.H"
|
|
#include "volFields.H"
|
|
#include "faceSet.H"
|
|
#include "cellSet.H"
|
|
#include "polyTopoChange.H"
|
|
#include "removeCells.H"
|
|
#include "edgeHashes.H"
|
|
#include "syncTools.H"
|
|
#include "ReadFields.H"
|
|
#include "mappedWallPolyPatch.H"
|
|
#include "fvMeshTools.H"
|
|
#include "zeroGradientFvPatchFields.H"
|
|
#include "processorMeshes.H"
|
|
|
|
using namespace Foam;
|
|
|
|
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
|
|
|
// Prepend prefix to selected patches.
|
|
void renamePatches
|
|
(
|
|
fvMesh& mesh,
|
|
const word& prefix,
|
|
const labelList& patchesToRename
|
|
)
|
|
{
|
|
polyBoundaryMesh& polyPatches =
|
|
const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
|
|
forAll(patchesToRename, i)
|
|
{
|
|
label patchi = patchesToRename[i];
|
|
polyPatch& pp = polyPatches[patchi];
|
|
|
|
if (isA<coupledPolyPatch>(pp))
|
|
{
|
|
WarningInFunction
|
|
<< "Encountered coupled patch " << pp.name()
|
|
<< ". Will only rename the patch itself,"
|
|
<< " not any referred patches."
|
|
<< " This might have to be done by hand."
|
|
<< endl;
|
|
}
|
|
|
|
pp.name() = prefix + '_' + pp.name();
|
|
}
|
|
// Recalculate any demand driven data (e.g. group to name lookup)
|
|
polyPatches.updateMesh();
|
|
}
|
|
|
|
|
|
template<class GeoField>
|
|
void subsetVolFields
|
|
(
|
|
const fvMesh& mesh,
|
|
const fvMesh& subMesh,
|
|
const labelList& cellMap,
|
|
const labelList& faceMap,
|
|
const labelHashSet& addedPatches
|
|
)
|
|
{
|
|
const labelList patchMap(identity(mesh.boundaryMesh().size()));
|
|
|
|
HashTable<const GeoField*> fields
|
|
(
|
|
mesh.objectRegistry::lookupClass<GeoField>()
|
|
);
|
|
forAllConstIters(fields, iter)
|
|
{
|
|
const GeoField& fld = *iter.val();
|
|
|
|
Info<< "Mapping field " << fld.name() << endl;
|
|
|
|
tmp<GeoField> tSubFld
|
|
(
|
|
fvMeshSubset::interpolate
|
|
(
|
|
fld,
|
|
subMesh,
|
|
patchMap,
|
|
cellMap,
|
|
faceMap
|
|
)
|
|
);
|
|
|
|
// Hack: set value to 0 for introduced patches (since don't
|
|
// get initialised.
|
|
forAll(tSubFld().boundaryField(), patchi)
|
|
{
|
|
if (addedPatches.found(patchi))
|
|
{
|
|
tSubFld.ref().boundaryFieldRef()[patchi] ==
|
|
typename GeoField::value_type(Zero);
|
|
}
|
|
}
|
|
|
|
// Store on subMesh
|
|
GeoField* subFld = tSubFld.ptr();
|
|
subFld->rename(fld.name());
|
|
subFld->writeOpt(IOobject::AUTO_WRITE);
|
|
subFld->store();
|
|
}
|
|
}
|
|
|
|
|
|
template<class GeoField>
|
|
void subsetSurfaceFields
|
|
(
|
|
const fvMesh& mesh,
|
|
const fvMesh& subMesh,
|
|
const labelList& cellMap,
|
|
const labelList& faceMap,
|
|
const labelHashSet& addedPatches
|
|
)
|
|
{
|
|
const labelList patchMap(identity(mesh.boundaryMesh().size()));
|
|
|
|
HashTable<const GeoField*> fields
|
|
(
|
|
mesh.objectRegistry::lookupClass<GeoField>()
|
|
);
|
|
forAllConstIters(fields, iter)
|
|
{
|
|
const GeoField& fld = *iter.val();
|
|
|
|
Info<< "Mapping field " << fld.name() << endl;
|
|
|
|
tmp<GeoField> tSubFld
|
|
(
|
|
fvMeshSubset::interpolate
|
|
(
|
|
fld,
|
|
subMesh,
|
|
patchMap,
|
|
cellMap,
|
|
faceMap
|
|
)
|
|
);
|
|
|
|
// Hack: set value to 0 for introduced patches (since don't
|
|
// get initialised.
|
|
forAll(tSubFld().boundaryField(), patchi)
|
|
{
|
|
if (addedPatches.found(patchi))
|
|
{
|
|
tSubFld.ref().boundaryFieldRef()[patchi] ==
|
|
typename GeoField::value_type(Zero);
|
|
}
|
|
}
|
|
|
|
// Store on subMesh
|
|
GeoField* subFld = tSubFld.ptr();
|
|
subFld->rename(fld.name());
|
|
subFld->writeOpt(IOobject::AUTO_WRITE);
|
|
subFld->store();
|
|
}
|
|
}
|
|
|
|
// Select all cells not in the region
|
|
labelList getNonRegionCells(const labelList& cellRegion, const label regionI)
|
|
{
|
|
DynamicList<label> nonRegionCells(cellRegion.size());
|
|
forAll(cellRegion, celli)
|
|
{
|
|
if (cellRegion[celli] != regionI)
|
|
{
|
|
nonRegionCells.append(celli);
|
|
}
|
|
}
|
|
return nonRegionCells.shrink();
|
|
}
|
|
|
|
|
|
void addToInterface
|
|
(
|
|
const polyMesh& mesh,
|
|
const label zoneID,
|
|
const label ownRegion,
|
|
const label neiRegion,
|
|
EdgeMap<Map<label>>& regionsToSize
|
|
)
|
|
{
|
|
edge interface
|
|
(
|
|
min(ownRegion, neiRegion),
|
|
max(ownRegion, neiRegion)
|
|
);
|
|
|
|
auto iter = regionsToSize.find(interface);
|
|
|
|
if (iter.found())
|
|
{
|
|
// Check if zone present
|
|
auto zoneIter = iter().find(zoneID);
|
|
if (zoneIter.found())
|
|
{
|
|
// Found zone. Increment count.
|
|
++(*zoneIter);
|
|
}
|
|
else
|
|
{
|
|
// New or no zone. Insert with count 1.
|
|
iter().insert(zoneID, 1);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Create new interface of size 1.
|
|
Map<label> zoneToSize;
|
|
zoneToSize.insert(zoneID, 1);
|
|
regionsToSize.insert(interface, zoneToSize);
|
|
}
|
|
}
|
|
|
|
|
|
// Get region-region interface name and sizes.
|
|
// Returns interfaces as straight list for looping in identical order.
|
|
void getInterfaceSizes
|
|
(
|
|
const polyMesh& mesh,
|
|
const bool useFaceZones,
|
|
const labelList& cellRegion,
|
|
const wordList& regionNames,
|
|
|
|
edgeList& interfaces,
|
|
List<Pair<word>>& interfaceNames,
|
|
labelList& interfaceSizes,
|
|
labelList& faceToInterface
|
|
)
|
|
{
|
|
// From region-region to faceZone (or -1) to number of faces.
|
|
|
|
EdgeMap<Map<label>> regionsToSize;
|
|
|
|
|
|
// Internal faces
|
|
// ~~~~~~~~~~~~~~
|
|
|
|
forAll(mesh.faceNeighbour(), facei)
|
|
{
|
|
label ownRegion = cellRegion[mesh.faceOwner()[facei]];
|
|
label neiRegion = cellRegion[mesh.faceNeighbour()[facei]];
|
|
|
|
if (ownRegion != neiRegion)
|
|
{
|
|
addToInterface
|
|
(
|
|
mesh,
|
|
(useFaceZones ? mesh.faceZones().whichZone(facei) : -1),
|
|
ownRegion,
|
|
neiRegion,
|
|
regionsToSize
|
|
);
|
|
}
|
|
}
|
|
|
|
// Boundary faces
|
|
// ~~~~~~~~~~~~~~
|
|
|
|
// Neighbour cellRegion.
|
|
labelList coupledRegion(mesh.nBoundaryFaces());
|
|
|
|
forAll(coupledRegion, i)
|
|
{
|
|
label celli = mesh.faceOwner()[i+mesh.nInternalFaces()];
|
|
coupledRegion[i] = cellRegion[celli];
|
|
}
|
|
syncTools::swapBoundaryFaceList(mesh, coupledRegion);
|
|
|
|
forAll(coupledRegion, i)
|
|
{
|
|
label facei = i+mesh.nInternalFaces();
|
|
label ownRegion = cellRegion[mesh.faceOwner()[facei]];
|
|
label neiRegion = coupledRegion[i];
|
|
|
|
if (ownRegion != neiRegion)
|
|
{
|
|
addToInterface
|
|
(
|
|
mesh,
|
|
(useFaceZones ? mesh.faceZones().whichZone(facei) : -1),
|
|
ownRegion,
|
|
neiRegion,
|
|
regionsToSize
|
|
);
|
|
}
|
|
}
|
|
|
|
|
|
if (Pstream::parRun())
|
|
{
|
|
if (Pstream::master())
|
|
{
|
|
// Receive and add to my sizes
|
|
for (const int slave : Pstream::subProcs())
|
|
{
|
|
IPstream fromSlave(Pstream::commsTypes::blocking, slave);
|
|
|
|
EdgeMap<Map<label>> slaveSizes(fromSlave);
|
|
|
|
forAllConstIters(slaveSizes, slaveIter)
|
|
{
|
|
const Map<label>& slaveInfo = *slaveIter;
|
|
|
|
auto masterIter = regionsToSize.find(slaveIter.key());
|
|
|
|
if (masterIter.found())
|
|
{
|
|
// Same inter-region
|
|
Map<label>& masterInfo = *masterIter;
|
|
|
|
forAllConstIters(slaveInfo, iter)
|
|
{
|
|
const label zoneID = iter.key();
|
|
const label slaveSize = iter.val();
|
|
|
|
auto zoneIter = masterInfo.find(zoneID);
|
|
if (zoneIter.found())
|
|
{
|
|
*zoneIter += slaveSize;
|
|
}
|
|
else
|
|
{
|
|
masterInfo.insert(zoneID, slaveSize);
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
regionsToSize.insert(slaveIter.key(), slaveInfo);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Send to master
|
|
{
|
|
OPstream toMaster
|
|
(
|
|
Pstream::commsTypes::blocking,
|
|
Pstream::masterNo()
|
|
);
|
|
toMaster << regionsToSize;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Rework
|
|
|
|
Pstream::scatter(regionsToSize);
|
|
|
|
|
|
|
|
// Now we have the global sizes of all inter-regions.
|
|
// Invert this on master and distribute.
|
|
label nInterfaces = 0;
|
|
forAllConstIters(regionsToSize, iter)
|
|
{
|
|
const Map<label>& info = iter.val();
|
|
nInterfaces += info.size();
|
|
}
|
|
|
|
interfaces.setSize(nInterfaces);
|
|
interfaceNames.setSize(nInterfaces);
|
|
interfaceSizes.setSize(nInterfaces);
|
|
EdgeMap<Map<label>> regionsToInterface(nInterfaces);
|
|
|
|
nInterfaces = 0;
|
|
forAllConstIters(regionsToSize, iter)
|
|
{
|
|
const edge& e = iter.key();
|
|
const Map<label>& info = iter.val();
|
|
|
|
const word& name0 = regionNames[e[0]];
|
|
const word& name1 = regionNames[e[1]];
|
|
|
|
forAllConstIters(info, infoIter)
|
|
{
|
|
interfaces[nInterfaces] = iter.key();
|
|
label zoneID = infoIter.key();
|
|
if (zoneID == -1)
|
|
{
|
|
interfaceNames[nInterfaces] = Pair<word>
|
|
(
|
|
name0 + "_to_" + name1,
|
|
name1 + "_to_" + name0
|
|
);
|
|
}
|
|
else
|
|
{
|
|
const word& zoneName = mesh.faceZones()[zoneID].name();
|
|
interfaceNames[nInterfaces] = Pair<word>
|
|
(
|
|
zoneName + "_" + name0 + "_to_" + name1,
|
|
zoneName + "_" + name1 + "_to_" + name0
|
|
);
|
|
}
|
|
interfaceSizes[nInterfaces] = infoIter();
|
|
|
|
if (regionsToInterface.found(e))
|
|
{
|
|
regionsToInterface[e].insert(zoneID, nInterfaces);
|
|
}
|
|
else
|
|
{
|
|
Map<label> zoneAndInterface;
|
|
zoneAndInterface.insert(zoneID, nInterfaces);
|
|
regionsToInterface.insert(e, zoneAndInterface);
|
|
}
|
|
nInterfaces++;
|
|
}
|
|
}
|
|
|
|
|
|
// Now all processor have consistent interface information
|
|
|
|
Pstream::scatter(interfaces);
|
|
Pstream::scatter(interfaceNames);
|
|
Pstream::scatter(interfaceSizes);
|
|
Pstream::scatter(regionsToInterface);
|
|
|
|
// Mark all inter-region faces.
|
|
faceToInterface.setSize(mesh.nFaces(), -1);
|
|
|
|
forAll(mesh.faceNeighbour(), facei)
|
|
{
|
|
label ownRegion = cellRegion[mesh.faceOwner()[facei]];
|
|
label neiRegion = cellRegion[mesh.faceNeighbour()[facei]];
|
|
|
|
if (ownRegion != neiRegion)
|
|
{
|
|
label zoneID = -1;
|
|
if (useFaceZones)
|
|
{
|
|
zoneID = mesh.faceZones().whichZone(facei);
|
|
}
|
|
|
|
edge interface
|
|
(
|
|
min(ownRegion, neiRegion),
|
|
max(ownRegion, neiRegion)
|
|
);
|
|
|
|
faceToInterface[facei] = regionsToInterface[interface][zoneID];
|
|
}
|
|
}
|
|
forAll(coupledRegion, i)
|
|
{
|
|
label facei = i+mesh.nInternalFaces();
|
|
label ownRegion = cellRegion[mesh.faceOwner()[facei]];
|
|
label neiRegion = coupledRegion[i];
|
|
|
|
if (ownRegion != neiRegion)
|
|
{
|
|
label zoneID = -1;
|
|
if (useFaceZones)
|
|
{
|
|
zoneID = mesh.faceZones().whichZone(facei);
|
|
}
|
|
|
|
edge interface
|
|
(
|
|
min(ownRegion, neiRegion),
|
|
max(ownRegion, neiRegion)
|
|
);
|
|
|
|
faceToInterface[facei] = regionsToInterface[interface][zoneID];
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Create mesh for region.
|
|
autoPtr<mapPolyMesh> createRegionMesh
|
|
(
|
|
const fvMesh& mesh,
|
|
// Region info
|
|
const labelList& cellRegion,
|
|
const label regionI,
|
|
const word& regionName,
|
|
// Interface info
|
|
const labelList& interfacePatches,
|
|
const labelList& faceToInterface,
|
|
|
|
autoPtr<fvMesh>& newMesh
|
|
)
|
|
{
|
|
// Create dummy system/fv*
|
|
fvMeshTools::createDummyFvMeshFiles(mesh, regionName, true);
|
|
|
|
// Neighbour cellRegion.
|
|
labelList coupledRegion(mesh.nBoundaryFaces());
|
|
|
|
forAll(coupledRegion, i)
|
|
{
|
|
label celli = mesh.faceOwner()[i+mesh.nInternalFaces()];
|
|
coupledRegion[i] = cellRegion[celli];
|
|
}
|
|
syncTools::swapBoundaryFaceList(mesh, coupledRegion);
|
|
|
|
|
|
// Topology change container. Start off from existing mesh.
|
|
polyTopoChange meshMod(mesh);
|
|
|
|
// Cell remover engine
|
|
removeCells cellRemover(mesh);
|
|
|
|
// Select all but region cells
|
|
labelList cellsToRemove(getNonRegionCells(cellRegion, regionI));
|
|
|
|
// Find out which faces will get exposed. Note that this
|
|
// gets faces in mesh face order. So both regions will get same
|
|
// face in same order (important!)
|
|
labelList exposedFaces = cellRemover.getExposedFaces(cellsToRemove);
|
|
|
|
labelList exposedPatchIDs(exposedFaces.size());
|
|
forAll(exposedFaces, i)
|
|
{
|
|
label facei = exposedFaces[i];
|
|
label interfacei = faceToInterface[facei];
|
|
|
|
label ownRegion = cellRegion[mesh.faceOwner()[facei]];
|
|
label neiRegion = -1;
|
|
|
|
if (mesh.isInternalFace(facei))
|
|
{
|
|
neiRegion = cellRegion[mesh.faceNeighbour()[facei]];
|
|
}
|
|
else
|
|
{
|
|
neiRegion = coupledRegion[facei-mesh.nInternalFaces()];
|
|
}
|
|
|
|
|
|
// Check which side is being kept - determines which of the two
|
|
// patches will be used.
|
|
|
|
label otherRegion = -1;
|
|
|
|
if (ownRegion == regionI && neiRegion != regionI)
|
|
{
|
|
otherRegion = neiRegion;
|
|
}
|
|
else if (ownRegion != regionI && neiRegion == regionI)
|
|
{
|
|
otherRegion = ownRegion;
|
|
}
|
|
else
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Exposed face:" << facei
|
|
<< " fc:" << mesh.faceCentres()[facei]
|
|
<< " has owner region " << ownRegion
|
|
<< " and neighbour region " << neiRegion
|
|
<< " when handling region:" << regionI
|
|
<< exit(FatalError);
|
|
}
|
|
|
|
// Find the patch.
|
|
if (regionI < otherRegion)
|
|
{
|
|
exposedPatchIDs[i] = interfacePatches[interfacei];
|
|
}
|
|
else
|
|
{
|
|
exposedPatchIDs[i] = interfacePatches[interfacei]+1;
|
|
}
|
|
}
|
|
|
|
// Remove faces
|
|
cellRemover.setRefinement
|
|
(
|
|
cellsToRemove,
|
|
exposedFaces,
|
|
exposedPatchIDs,
|
|
meshMod
|
|
);
|
|
|
|
autoPtr<mapPolyMesh> map = meshMod.makeMesh
|
|
(
|
|
newMesh,
|
|
IOobject
|
|
(
|
|
regionName,
|
|
mesh.time().timeName(),
|
|
mesh.time(),
|
|
IOobject::READ_IF_PRESENT, // read fv* if present
|
|
IOobject::AUTO_WRITE
|
|
),
|
|
mesh
|
|
);
|
|
|
|
return map;
|
|
}
|
|
|
|
|
|
void createAndWriteRegion
|
|
(
|
|
const fvMesh& mesh,
|
|
const labelList& cellRegion,
|
|
const wordList& regionNames,
|
|
const bool prefixRegion,
|
|
const labelList& faceToInterface,
|
|
const labelList& interfacePatches,
|
|
const label regionI,
|
|
const word& newMeshInstance
|
|
)
|
|
{
|
|
Info<< "Creating mesh for region " << regionI
|
|
<< ' ' << regionNames[regionI] << endl;
|
|
|
|
autoPtr<fvMesh> newMesh;
|
|
autoPtr<mapPolyMesh> map = createRegionMesh
|
|
(
|
|
mesh,
|
|
cellRegion,
|
|
regionI,
|
|
regionNames[regionI],
|
|
interfacePatches,
|
|
faceToInterface,
|
|
newMesh
|
|
);
|
|
|
|
|
|
// Make map of all added patches
|
|
labelHashSet addedPatches(2*interfacePatches.size());
|
|
forAll(interfacePatches, interfacei)
|
|
{
|
|
addedPatches.insert(interfacePatches[interfacei]);
|
|
addedPatches.insert(interfacePatches[interfacei]+1);
|
|
}
|
|
|
|
|
|
Info<< "Mapping fields" << endl;
|
|
|
|
// Map existing fields
|
|
newMesh().updateMesh(map());
|
|
|
|
// Add subsetted fields
|
|
subsetVolFields<volScalarField>
|
|
(
|
|
mesh,
|
|
newMesh(),
|
|
map().cellMap(),
|
|
map().faceMap(),
|
|
addedPatches
|
|
);
|
|
subsetVolFields<volVectorField>
|
|
(
|
|
mesh,
|
|
newMesh(),
|
|
map().cellMap(),
|
|
map().faceMap(),
|
|
addedPatches
|
|
);
|
|
subsetVolFields<volSphericalTensorField>
|
|
(
|
|
mesh,
|
|
newMesh(),
|
|
map().cellMap(),
|
|
map().faceMap(),
|
|
addedPatches
|
|
);
|
|
subsetVolFields<volSymmTensorField>
|
|
(
|
|
mesh,
|
|
newMesh(),
|
|
map().cellMap(),
|
|
map().faceMap(),
|
|
addedPatches
|
|
);
|
|
subsetVolFields<volTensorField>
|
|
(
|
|
mesh,
|
|
newMesh(),
|
|
map().cellMap(),
|
|
map().faceMap(),
|
|
addedPatches
|
|
);
|
|
|
|
subsetSurfaceFields<surfaceScalarField>
|
|
(
|
|
mesh,
|
|
newMesh(),
|
|
map().cellMap(),
|
|
map().faceMap(),
|
|
addedPatches
|
|
);
|
|
subsetSurfaceFields<surfaceVectorField>
|
|
(
|
|
mesh,
|
|
newMesh(),
|
|
map().cellMap(),
|
|
map().faceMap(),
|
|
addedPatches
|
|
);
|
|
subsetSurfaceFields<surfaceSphericalTensorField>
|
|
(
|
|
mesh,
|
|
newMesh(),
|
|
map().cellMap(),
|
|
map().faceMap(),
|
|
addedPatches
|
|
);
|
|
subsetSurfaceFields<surfaceSymmTensorField>
|
|
(
|
|
mesh,
|
|
newMesh(),
|
|
map().cellMap(),
|
|
map().faceMap(),
|
|
addedPatches
|
|
);
|
|
subsetSurfaceFields<surfaceTensorField>
|
|
(
|
|
mesh,
|
|
newMesh(),
|
|
map().cellMap(),
|
|
map().faceMap(),
|
|
addedPatches
|
|
);
|
|
|
|
|
|
const polyBoundaryMesh& newPatches = newMesh().boundaryMesh();
|
|
newPatches.checkParallelSync(true);
|
|
|
|
// Delete empty patches
|
|
// ~~~~~~~~~~~~~~~~~~~~
|
|
|
|
// Create reordering list to move patches-to-be-deleted to end
|
|
labelList oldToNew(newPatches.size(), -1);
|
|
DynamicList<label> sharedPatches(newPatches.size());
|
|
label newI = 0;
|
|
|
|
Info<< "Deleting empty patches" << endl;
|
|
|
|
// Assumes all non-proc boundaries are on all processors!
|
|
forAll(newPatches, patchi)
|
|
{
|
|
const polyPatch& pp = newPatches[patchi];
|
|
|
|
if (!isA<processorPolyPatch>(pp))
|
|
{
|
|
if (returnReduce(pp.size(), sumOp<label>()) > 0)
|
|
{
|
|
oldToNew[patchi] = newI;
|
|
if (!addedPatches.found(patchi))
|
|
{
|
|
sharedPatches.append(newI);
|
|
}
|
|
newI++;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Same for processor patches (but need no reduction)
|
|
forAll(newPatches, patchi)
|
|
{
|
|
const polyPatch& pp = newPatches[patchi];
|
|
|
|
if (isA<processorPolyPatch>(pp) && pp.size())
|
|
{
|
|
oldToNew[patchi] = newI++;
|
|
}
|
|
}
|
|
|
|
const label nNewPatches = newI;
|
|
|
|
// Move all deleteable patches to the end
|
|
forAll(oldToNew, patchi)
|
|
{
|
|
if (oldToNew[patchi] == -1)
|
|
{
|
|
oldToNew[patchi] = newI++;
|
|
}
|
|
}
|
|
|
|
//reorderPatches(newMesh(), oldToNew, nNewPatches);
|
|
fvMeshTools::reorderPatches(newMesh(), oldToNew, nNewPatches, true);
|
|
|
|
// Rename shared patches with region name
|
|
if (prefixRegion)
|
|
{
|
|
Info<< "Prefixing patches with region name" << endl;
|
|
|
|
renamePatches(newMesh(), regionNames[regionI], sharedPatches);
|
|
}
|
|
|
|
|
|
Info<< "Writing new mesh" << endl;
|
|
|
|
newMesh().setInstance(newMeshInstance);
|
|
newMesh().write();
|
|
topoSet::removeFiles(newMesh());
|
|
processorMeshes::removeFiles(newMesh());
|
|
|
|
// Write addressing files like decomposePar
|
|
Info<< "Writing addressing to base mesh" << endl;
|
|
|
|
labelIOList pointProcAddressing
|
|
(
|
|
IOobject
|
|
(
|
|
"pointRegionAddressing",
|
|
newMesh().facesInstance(),
|
|
newMesh().meshSubDir,
|
|
newMesh(),
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
),
|
|
map().pointMap()
|
|
);
|
|
Info<< "Writing map " << pointProcAddressing.name()
|
|
<< " from region" << regionI
|
|
<< " points back to base mesh." << endl;
|
|
pointProcAddressing.write();
|
|
|
|
labelIOList faceProcAddressing
|
|
(
|
|
IOobject
|
|
(
|
|
"faceRegionAddressing",
|
|
newMesh().facesInstance(),
|
|
newMesh().meshSubDir,
|
|
newMesh(),
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
),
|
|
newMesh().nFaces()
|
|
);
|
|
forAll(faceProcAddressing, facei)
|
|
{
|
|
// face + turning index. (see decomposePar)
|
|
// Is the face pointing in the same direction?
|
|
label oldFacei = map().faceMap()[facei];
|
|
|
|
if
|
|
(
|
|
map().cellMap()[newMesh().faceOwner()[facei]]
|
|
== mesh.faceOwner()[oldFacei]
|
|
)
|
|
{
|
|
faceProcAddressing[facei] = oldFacei+1;
|
|
}
|
|
else
|
|
{
|
|
faceProcAddressing[facei] = -(oldFacei+1);
|
|
}
|
|
}
|
|
Info<< "Writing map " << faceProcAddressing.name()
|
|
<< " from region" << regionI
|
|
<< " faces back to base mesh." << endl;
|
|
faceProcAddressing.write();
|
|
|
|
labelIOList cellProcAddressing
|
|
(
|
|
IOobject
|
|
(
|
|
"cellRegionAddressing",
|
|
newMesh().facesInstance(),
|
|
newMesh().meshSubDir,
|
|
newMesh(),
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
),
|
|
map().cellMap()
|
|
);
|
|
Info<< "Writing map " <<cellProcAddressing.name()
|
|
<< " from region" << regionI
|
|
<< " cells back to base mesh." << endl;
|
|
cellProcAddressing.write();
|
|
|
|
labelIOList boundaryProcAddressing
|
|
(
|
|
IOobject
|
|
(
|
|
"boundaryRegionAddressing",
|
|
newMesh().facesInstance(),
|
|
newMesh().meshSubDir,
|
|
newMesh(),
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
),
|
|
labelList(nNewPatches, -1)
|
|
);
|
|
forAll(oldToNew, i)
|
|
{
|
|
if (!addedPatches.found(i))
|
|
{
|
|
label newI = oldToNew[i];
|
|
if (newI >= 0 && newI < nNewPatches)
|
|
{
|
|
boundaryProcAddressing[oldToNew[i]] = i;
|
|
}
|
|
}
|
|
}
|
|
Info<< "Writing map " << boundaryProcAddressing.name()
|
|
<< " from region" << regionI
|
|
<< " boundary back to base mesh." << endl;
|
|
boundaryProcAddressing.write();
|
|
}
|
|
|
|
|
|
// Create for every region-region interface with non-zero size two patches.
|
|
// First one is for minimumregion to maximumregion.
|
|
// Note that patches get created in same order on all processors (if parallel)
|
|
// since looping over synchronised 'interfaces'.
|
|
labelList addRegionPatches
|
|
(
|
|
fvMesh& mesh,
|
|
const wordList& regionNames,
|
|
const edgeList& interfaces,
|
|
const List<Pair<word>>& interfaceNames
|
|
)
|
|
{
|
|
Info<< nl << "Adding patches" << nl << endl;
|
|
|
|
labelList interfacePatches(interfaces.size());
|
|
|
|
forAll(interfaces, interI)
|
|
{
|
|
const edge& e = interfaces[interI];
|
|
const Pair<word>& names = interfaceNames[interI];
|
|
|
|
//Info<< "For interface " << interI
|
|
// << " between regions " << e
|
|
// << " trying to add patches " << names << endl;
|
|
|
|
|
|
mappedWallPolyPatch patch1
|
|
(
|
|
names[0],
|
|
0, // overridden
|
|
0, // overridden
|
|
0, // overridden
|
|
regionNames[e[1]], // sampleRegion
|
|
mappedPatchBase::NEARESTPATCHFACE,
|
|
names[1], // samplePatch
|
|
point::zero, // offset
|
|
mesh.boundaryMesh()
|
|
);
|
|
|
|
interfacePatches[interI] = fvMeshTools::addPatch
|
|
(
|
|
mesh,
|
|
patch1,
|
|
dictionary(), //optional per field value
|
|
calculatedFvPatchField<scalar>::typeName,
|
|
true //validBoundary
|
|
);
|
|
|
|
mappedWallPolyPatch patch2
|
|
(
|
|
names[1],
|
|
0,
|
|
0,
|
|
0,
|
|
regionNames[e[0]], // sampleRegion
|
|
mappedPatchBase::NEARESTPATCHFACE,
|
|
names[0],
|
|
point::zero, // offset
|
|
mesh.boundaryMesh()
|
|
);
|
|
fvMeshTools::addPatch
|
|
(
|
|
mesh,
|
|
patch2,
|
|
dictionary(), //optional per field value
|
|
calculatedFvPatchField<scalar>::typeName,
|
|
true //validBoundary
|
|
);
|
|
|
|
Info<< "For interface between region " << regionNames[e[0]]
|
|
<< " and " << regionNames[e[1]] << " added patches" << endl
|
|
<< " " << interfacePatches[interI]
|
|
<< "\t" << mesh.boundaryMesh()[interfacePatches[interI]].name()
|
|
<< endl
|
|
<< " " << interfacePatches[interI]+1
|
|
<< "\t" << mesh.boundaryMesh()[interfacePatches[interI]+1].name()
|
|
<< endl;
|
|
}
|
|
return interfacePatches;
|
|
}
|
|
|
|
|
|
// Find region that covers most of cell zone
|
|
label findCorrespondingRegion
|
|
(
|
|
const labelList& existingZoneID, // per cell the (unique) zoneID
|
|
const labelList& cellRegion,
|
|
const label nCellRegions,
|
|
const label zoneI,
|
|
const label minOverlapSize
|
|
)
|
|
{
|
|
// Per region the number of cells in zoneI
|
|
labelList cellsInZone(nCellRegions, Zero);
|
|
|
|
forAll(cellRegion, celli)
|
|
{
|
|
if (existingZoneID[celli] == zoneI)
|
|
{
|
|
cellsInZone[cellRegion[celli]]++;
|
|
}
|
|
}
|
|
|
|
Pstream::listCombineGather(cellsInZone, plusEqOp<label>());
|
|
Pstream::listCombineScatter(cellsInZone);
|
|
|
|
// Pick region with largest overlap of zoneI
|
|
label regionI = findMax(cellsInZone);
|
|
|
|
|
|
if (cellsInZone[regionI] < minOverlapSize)
|
|
{
|
|
// Region covers too little of zone. Not good enough.
|
|
regionI = -1;
|
|
}
|
|
else
|
|
{
|
|
// Check that region contains no cells that aren't in cellZone.
|
|
forAll(cellRegion, celli)
|
|
{
|
|
if (cellRegion[celli] == regionI && existingZoneID[celli] != zoneI)
|
|
{
|
|
// celli in regionI but not in zoneI
|
|
regionI = -1;
|
|
break;
|
|
}
|
|
}
|
|
// If one in error, all should be in error. Note that branch gets taken
|
|
// on all procs.
|
|
reduce(regionI, minOp<label>());
|
|
}
|
|
|
|
return regionI;
|
|
}
|
|
|
|
|
|
void getClusterID
|
|
(
|
|
const polyMesh& mesh,
|
|
const cellZoneMesh& cellZones,
|
|
const wordList& clusterNames,
|
|
const labelListList& clusterToZones,
|
|
labelList& clusterID,
|
|
labelList& neiClusterID
|
|
)
|
|
{
|
|
// Existing zoneID
|
|
clusterID.setSize(mesh.nCells());
|
|
clusterID = -1;
|
|
|
|
forAll(clusterToZones, clusterI)
|
|
{
|
|
for (const label zoneI : clusterToZones[clusterI])
|
|
{
|
|
const cellZone& cz = cellZones[zoneI];
|
|
|
|
forAll(cz, i)
|
|
{
|
|
label celli = cz[i];
|
|
if (clusterID[celli] == -1)
|
|
{
|
|
clusterID[celli] = clusterI;
|
|
}
|
|
else
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Cell " << celli << " with cell centre "
|
|
<< mesh.cellCentres()[celli]
|
|
<< " is multiple zones. This is not allowed." << endl
|
|
<< "It is in zone " << clusterNames[clusterID[celli]]
|
|
<< " and in zone " << clusterNames[clusterI]
|
|
<< exit(FatalError);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Neighbour zoneID.
|
|
syncTools::swapBoundaryCellList(mesh, clusterID, neiClusterID);
|
|
}
|
|
|
|
|
|
word makeRegionName
|
|
(
|
|
const cellZoneMesh& czs,
|
|
const label regioni,
|
|
const labelList& zoneIDs
|
|
)
|
|
{
|
|
// Synthesise region name. Equals the zone name if cluster consist of only
|
|
// one zone
|
|
|
|
if (zoneIDs.empty())
|
|
{
|
|
return word("domain") + Foam::name(regioni);
|
|
}
|
|
else
|
|
{
|
|
// Synthesize name from appended cellZone names
|
|
word regionName(czs[zoneIDs[0]].name());
|
|
for (label i = 1; i < zoneIDs.size(); i++)
|
|
{
|
|
regionName += "_" + czs[zoneIDs[i]].name();
|
|
}
|
|
return regionName;
|
|
}
|
|
}
|
|
|
|
|
|
void makeClusters
|
|
(
|
|
const List<wordRes>& zoneClusters,
|
|
const cellZoneMesh& cellZones,
|
|
wordList& clusterNames,
|
|
labelListList& clusterToZones,
|
|
labelList& zoneToCluster
|
|
)
|
|
{
|
|
// Check if there are clustering for zones. If none every zone goes into
|
|
// its own cluster.
|
|
|
|
clusterNames.clear();
|
|
clusterToZones.clear();
|
|
zoneToCluster.setSize(cellZones.size());
|
|
zoneToCluster = -1;
|
|
|
|
if (zoneClusters.size())
|
|
{
|
|
forAll(zoneClusters, clusteri)
|
|
{
|
|
const labelList zoneIDs(cellZones.indices(zoneClusters[clusteri]));
|
|
UIndirectList<label>(zoneToCluster, zoneIDs) = clusteri;
|
|
clusterNames.append(makeRegionName(cellZones, clusteri, zoneIDs));
|
|
clusterToZones.append(std::move(zoneIDs));
|
|
}
|
|
|
|
// Unclustered zone
|
|
forAll(zoneToCluster, zonei)
|
|
{
|
|
if (zoneToCluster[zonei] == -1)
|
|
{
|
|
clusterNames.append(cellZones[zonei].name());
|
|
clusterToZones.append(labelList(1, zonei));
|
|
zoneToCluster[zonei] = clusterToZones.size();
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
for (const auto& cellZone : cellZones)
|
|
{
|
|
const label nClusters = clusterToZones.size();
|
|
clusterNames.append(cellZone.name());
|
|
clusterToZones.append(labelList(1, cellZone.index()));
|
|
zoneToCluster[cellZone.index()] = nClusters;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void matchRegions
|
|
(
|
|
const bool sloppyCellZones,
|
|
const polyMesh& mesh,
|
|
|
|
const wordList& clusterNames,
|
|
const labelListList& clusterToZones,
|
|
const labelList& clusterID,
|
|
|
|
const label nCellRegions,
|
|
const labelList& cellRegion,
|
|
|
|
labelListList& regionToZones,
|
|
wordList& regionNames,
|
|
labelList& zoneToRegion
|
|
)
|
|
{
|
|
const cellZoneMesh& cellZones = mesh.cellZones();
|
|
|
|
regionToZones.setSize(nCellRegions);
|
|
regionToZones = labelList();
|
|
regionNames.setSize(nCellRegions);
|
|
regionNames = word::null;
|
|
zoneToRegion.setSize(cellZones.size(), -1);
|
|
|
|
|
|
// Sizes per cluster
|
|
labelList clusterSizes(clusterToZones.size(), Zero);
|
|
forAll(clusterToZones, clusterI)
|
|
{
|
|
for (const label zoneI : clusterToZones[clusterI])
|
|
{
|
|
clusterSizes[clusterI] += cellZones[zoneI].size();
|
|
}
|
|
reduce(clusterSizes[clusterI], sumOp<label>());
|
|
}
|
|
|
|
if (sloppyCellZones)
|
|
{
|
|
Info<< "Trying to match regions to existing cell zones;"
|
|
<< " region can be subset of cell zone." << nl << endl;
|
|
|
|
forAll(clusterToZones, clusterI)
|
|
{
|
|
label regionI = findCorrespondingRegion
|
|
(
|
|
clusterID,
|
|
cellRegion,
|
|
nCellRegions,
|
|
clusterI,
|
|
label(0.5*clusterSizes[clusterI]) // minimum overlap
|
|
);
|
|
|
|
if (regionI != -1)
|
|
{
|
|
Info<< "Sloppily matched region " << regionI
|
|
//<< " size " << regionSizes[regionI]
|
|
<< " to cluster " << clusterI
|
|
<< " size " << clusterSizes[clusterI]
|
|
<< endl;
|
|
UIndirectList<label>
|
|
(
|
|
zoneToRegion,
|
|
clusterToZones[clusterI]
|
|
) = regionI;
|
|
regionToZones[regionI] = clusterToZones[clusterI];
|
|
regionNames[regionI] = clusterNames[clusterI];
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
Info<< "Trying to match regions to existing cell zones." << nl << endl;
|
|
|
|
forAll(clusterToZones, clusterI)
|
|
{
|
|
label regionI = findCorrespondingRegion
|
|
(
|
|
clusterID,
|
|
cellRegion,
|
|
nCellRegions,
|
|
clusterI,
|
|
clusterSizes[clusterI] // require exact match
|
|
);
|
|
|
|
if (regionI != -1)
|
|
{
|
|
UIndirectList<label>
|
|
(
|
|
zoneToRegion,
|
|
clusterToZones[clusterI]
|
|
) = regionI;
|
|
regionToZones[regionI] = clusterToZones[clusterI];
|
|
regionNames[regionI] = clusterNames[clusterI];
|
|
}
|
|
}
|
|
}
|
|
// Allocate region names for unmatched regions.
|
|
forAll(regionNames, regionI)
|
|
{
|
|
regionNames[regionI] = makeRegionName
|
|
(
|
|
cellZones,
|
|
regionI,
|
|
regionToZones[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(dimless, Zero),
|
|
zeroGradientFvPatchScalarField::typeName
|
|
);
|
|
forAll(cellRegion, celli)
|
|
{
|
|
cellToRegion[celli] = cellRegion[celli];
|
|
}
|
|
cellToRegion.write();
|
|
|
|
Info<< "Writing region per cell as volScalarField to "
|
|
<< cellToRegion.objectPath() << nl << endl;
|
|
}
|
|
}
|
|
|
|
|
|
|
|
|
|
int main(int argc, char *argv[])
|
|
{
|
|
argList::addNote
|
|
(
|
|
"Split mesh into multiple regions (detected by walking across faces)"
|
|
);
|
|
#include "addRegionOption.H"
|
|
#include "addOverwriteOption.H"
|
|
argList::addBoolOption
|
|
(
|
|
"cellZones",
|
|
"Additionally split cellZones off into separate regions"
|
|
);
|
|
argList::addBoolOption
|
|
(
|
|
"cellZonesOnly",
|
|
"Use cellZones only to split mesh into regions; do not use walking"
|
|
);
|
|
argList::addOption
|
|
(
|
|
"cellZonesFileOnly",
|
|
"file",
|
|
"Like -cellZonesOnly, but use specified file"
|
|
);
|
|
argList::addOption
|
|
(
|
|
"combineZones",
|
|
"lists of zones",
|
|
"Combine zones in follow-on analysis"
|
|
);
|
|
argList::addOption
|
|
(
|
|
"blockedFaces",
|
|
"faceSet",
|
|
"Specify additional region boundaries that walking does not cross"
|
|
);
|
|
argList::addBoolOption
|
|
(
|
|
"makeCellZones",
|
|
"Place cells into cellZones instead of splitting mesh"
|
|
);
|
|
argList::addBoolOption
|
|
(
|
|
"largestOnly",
|
|
"Only write largest region"
|
|
);
|
|
argList::addOption
|
|
(
|
|
"insidePoint",
|
|
"point",
|
|
"Only write region containing point"
|
|
);
|
|
argList::addBoolOption
|
|
(
|
|
"detectOnly",
|
|
"Do not write mesh"
|
|
);
|
|
argList::addBoolOption
|
|
(
|
|
"sloppyCellZones",
|
|
"Try to match heuristically regions to existing cell zones"
|
|
);
|
|
argList::addBoolOption
|
|
(
|
|
"useFaceZones",
|
|
"Use faceZones to patch inter-region faces instead of single patch"
|
|
);
|
|
argList::addBoolOption
|
|
(
|
|
"prefixRegion",
|
|
"Prefix region name to all patches, not just coupling patches"
|
|
);
|
|
|
|
argList::noFunctionObjects(); // Never use function objects
|
|
|
|
#include "setRootCase.H"
|
|
#include "createTime.H"
|
|
#include "createNamedMesh.H"
|
|
|
|
const word oldInstance = mesh.pointsInstance();
|
|
|
|
word blockedFacesName;
|
|
if (args.readIfPresent("blockedFaces", blockedFacesName))
|
|
{
|
|
Info<< "Reading blocked internal faces from faceSet "
|
|
<< blockedFacesName << nl << endl;
|
|
}
|
|
|
|
const bool makeCellZones = args.found("makeCellZones");
|
|
const bool largestOnly = args.found("largestOnly");
|
|
const bool insidePoint = args.found("insidePoint");
|
|
const bool useCellZones = args.found("cellZones");
|
|
const bool useCellZonesOnly = args.found("cellZonesOnly");
|
|
const bool useCellZonesFile = args.found("cellZonesFileOnly");
|
|
const bool combineZones = args.found("combineZones");
|
|
const bool overwrite = args.found("overwrite");
|
|
const bool detectOnly = args.found("detectOnly");
|
|
const bool sloppyCellZones = args.found("sloppyCellZones");
|
|
const bool useFaceZones = args.found("useFaceZones");
|
|
const bool prefixRegion = args.found("prefixRegion");
|
|
|
|
|
|
if
|
|
(
|
|
(useCellZonesOnly || useCellZonesFile)
|
|
&& (useCellZones || blockedFacesName.size())
|
|
)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "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 (useFaceZones)
|
|
{
|
|
Info<< "Using current faceZones to divide inter-region interfaces"
|
|
<< " into multiple patches."
|
|
<< nl << endl;
|
|
}
|
|
else
|
|
{
|
|
Info<< "Creating single patch per inter-region interface."
|
|
<< nl << endl;
|
|
}
|
|
|
|
|
|
|
|
if (insidePoint && largestOnly)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "You cannot specify both -largestOnly"
|
|
<< " (keep region with most cells)"
|
|
<< " and -insidePoint (keep region containing point)"
|
|
<< exit(FatalError);
|
|
}
|
|
|
|
|
|
// Make sure cellZone names consistent across processors
|
|
mesh.cellZones().checkParallelSync(true);
|
|
|
|
List<wordRes> zoneClusters;
|
|
if (combineZones)
|
|
{
|
|
zoneClusters = args.get<List<wordRes>>("combineZones");
|
|
}
|
|
|
|
|
|
// Determine per cell the region it belongs to
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
// cellRegion is the labelList with the region per cell.
|
|
labelList cellRegion;
|
|
// Region to zone(s)
|
|
labelListList regionToZones;
|
|
// 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;
|
|
|
|
// Collect sets of zones into clusters. If no cluster is just an identity
|
|
// list (cluster 0 is cellZone 0 etc.)
|
|
wordList clusterNames;
|
|
labelListList clusterToZones;
|
|
labelList zoneToCluster;
|
|
makeClusters
|
|
(
|
|
zoneClusters,
|
|
mesh.cellZones(),
|
|
clusterNames,
|
|
clusterToZones,
|
|
zoneToCluster
|
|
);
|
|
|
|
// Existing clusterID
|
|
labelList clusterID(mesh.nCells(), -1);
|
|
// Neighbour clusterID.
|
|
labelList neiClusterID(mesh.nBoundaryFaces());
|
|
getClusterID
|
|
(
|
|
mesh,
|
|
mesh.cellZones(),
|
|
clusterNames,
|
|
clusterToZones,
|
|
clusterID,
|
|
neiClusterID
|
|
);
|
|
|
|
label unzonedCelli = clusterID.find(-1);
|
|
if (unzonedCelli != -1)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "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 = clusterID;
|
|
nCellRegions = gMax(cellRegion)+1;
|
|
zoneToRegion = zoneToCluster;
|
|
regionToZones = clusterToZones;
|
|
regionNames = clusterNames;
|
|
}
|
|
else if (useCellZonesFile)
|
|
{
|
|
const word zoneFile(args["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
|
|
);
|
|
|
|
wordList clusterNames;
|
|
labelListList clusterToZones;
|
|
labelList zoneToCluster;
|
|
makeClusters
|
|
(
|
|
zoneClusters,
|
|
newCellZones,
|
|
clusterNames,
|
|
clusterToZones,
|
|
zoneToCluster
|
|
);
|
|
|
|
|
|
// Existing clusterID
|
|
labelList clusterID(mesh.nCells(), -1);
|
|
// Neighbour clusterID.
|
|
labelList neiClusterID(mesh.nBoundaryFaces());
|
|
getClusterID
|
|
(
|
|
mesh,
|
|
newCellZones,
|
|
clusterNames,
|
|
clusterToZones,
|
|
clusterID,
|
|
neiClusterID
|
|
);
|
|
|
|
|
|
label unzonedCelli = clusterID.find(-1);
|
|
if (unzonedCelli != -1)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "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 = clusterID;
|
|
nCellRegions = gMax(cellRegion)+1;
|
|
zoneToRegion = zoneToCluster;
|
|
regionToZones = clusterToZones;
|
|
regionNames = clusterNames;
|
|
}
|
|
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);
|
|
|
|
for (const label facei : blockedFaceSet)
|
|
{
|
|
blockedFace[facei] = true;
|
|
}
|
|
}
|
|
|
|
// Collect sets of zones into clusters. If no cluster is just an
|
|
// identity list (cluster 0 is cellZone 0 etc.)
|
|
wordList clusterNames;
|
|
labelListList clusterToZones;
|
|
labelList zoneToCluster;
|
|
makeClusters
|
|
(
|
|
zoneClusters,
|
|
mesh.cellZones(),
|
|
clusterNames,
|
|
clusterToZones,
|
|
zoneToCluster
|
|
);
|
|
|
|
// Existing clusterID
|
|
labelList clusterID(mesh.nCells(), -1);
|
|
// Neighbour clusterID.
|
|
labelList neiClusterID(mesh.nBoundaryFaces());
|
|
getClusterID
|
|
(
|
|
mesh,
|
|
mesh.cellZones(),
|
|
clusterNames,
|
|
clusterToZones,
|
|
clusterID,
|
|
neiClusterID
|
|
);
|
|
|
|
|
|
// Imply from differing cellZones
|
|
if (useCellZones)
|
|
{
|
|
blockedFace.setSize(mesh.nFaces(), false);
|
|
|
|
for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
|
|
{
|
|
label ownCluster = clusterID[mesh.faceOwner()[facei]];
|
|
label neiCluster = clusterID[mesh.faceNeighbour()[facei]];
|
|
|
|
if (ownCluster != neiCluster)
|
|
{
|
|
blockedFace[facei] = true;
|
|
}
|
|
}
|
|
|
|
// Different cellZones on either side of processor patch.
|
|
forAll(neiClusterID, i)
|
|
{
|
|
label facei = i+mesh.nInternalFaces();
|
|
label ownCluster = clusterID[mesh.faceOwner()[facei]];
|
|
label neiCluster = neiClusterID[i];
|
|
|
|
if (ownCluster != neiCluster)
|
|
{
|
|
blockedFace[facei] = true;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Do a topological walk to determine regions
|
|
regionSplit regions(mesh, blockedFace);
|
|
nCellRegions = regions.nRegions();
|
|
cellRegion.transfer(regions);
|
|
|
|
// Match regions to zones
|
|
matchRegions
|
|
(
|
|
sloppyCellZones,
|
|
mesh,
|
|
|
|
// cluster-to-zone and cluster-to-cell addressing
|
|
clusterNames,
|
|
clusterToZones,
|
|
clusterID,
|
|
|
|
// cell-to-region addressing
|
|
nCellRegions,
|
|
cellRegion,
|
|
|
|
// matched zones
|
|
regionToZones,
|
|
regionNames,
|
|
zoneToRegion
|
|
);
|
|
|
|
// Override any default region names if single region selected
|
|
if (largestOnly || insidePoint)
|
|
{
|
|
forAll(regionToZones, regionI)
|
|
{
|
|
if (regionToZones[regionI].empty())
|
|
{
|
|
if (overwrite)
|
|
{
|
|
regionNames[regionI] = polyMesh::defaultRegion;
|
|
}
|
|
else if (insidePoint)
|
|
{
|
|
regionNames[regionI] = "insidePoint";
|
|
}
|
|
else if (largestOnly)
|
|
{
|
|
regionNames[regionI] = "largestOnly";
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
Info<< endl << "Number of regions:" << nCellRegions << nl << endl;
|
|
|
|
|
|
// Write decomposition to file
|
|
writeCellToRegion(mesh, cellRegion);
|
|
|
|
|
|
|
|
// Sizes per region
|
|
// ~~~~~~~~~~~~~~~~
|
|
|
|
labelList regionSizes(nCellRegions, Zero);
|
|
|
|
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\t" << regionSizes[regionI] << nl;
|
|
}
|
|
Info<< endl;
|
|
|
|
|
|
|
|
// Print region to zone
|
|
Info<< "Region\tZone\tName" << nl
|
|
<< "------\t----\t----" << endl;
|
|
forAll(regionToZones, regionI)
|
|
{
|
|
Info<< regionI << "\t\t" << flatOutput(regionToZones[regionI])
|
|
<< '\t'
|
|
<< regionNames[regionI] << nl;
|
|
}
|
|
Info<< endl;
|
|
|
|
|
|
|
|
// Since we're going to mess with patches and zones make sure all
|
|
// is synchronised
|
|
mesh.boundaryMesh().checkParallelSync(true);
|
|
mesh.faceZones().checkParallelSync(true);
|
|
|
|
|
|
// Interfaces
|
|
// ----------
|
|
// per interface:
|
|
// - the two regions on either side
|
|
// - the name
|
|
// - the (global) size
|
|
edgeList interfaces;
|
|
List<Pair<word>> interfaceNames;
|
|
labelList interfaceSizes;
|
|
// per face the interface
|
|
labelList faceToInterface;
|
|
|
|
getInterfaceSizes
|
|
(
|
|
mesh,
|
|
useFaceZones,
|
|
cellRegion,
|
|
regionNames,
|
|
|
|
interfaces,
|
|
interfaceNames,
|
|
interfaceSizes,
|
|
faceToInterface
|
|
);
|
|
|
|
Info<< "Sizes of interfaces between regions:" << nl << nl
|
|
<< "Interface\tRegion\tRegion\tFaces" << nl
|
|
<< "---------\t------\t------\t-----" << endl;
|
|
|
|
forAll(interfaces, interI)
|
|
{
|
|
const edge& e = interfaces[interI];
|
|
|
|
Info<< interI
|
|
<< "\t\t\t" << e[0] << "\t\t" << e[1]
|
|
<< "\t\t" << interfaceSizes[interI] << nl;
|
|
}
|
|
Info<< endl;
|
|
|
|
|
|
if (detectOnly)
|
|
{
|
|
Info<< "End\n" << endl;
|
|
|
|
return 0;
|
|
}
|
|
|
|
|
|
// Read objects in time directory
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
IOobjectList objects(mesh, runTime.timeName());
|
|
|
|
// Read vol fields.
|
|
|
|
PtrList<volScalarField> vsFlds;
|
|
ReadFields(mesh, objects, vsFlds);
|
|
|
|
PtrList<volVectorField> vvFlds;
|
|
ReadFields(mesh, objects, vvFlds);
|
|
|
|
PtrList<volSphericalTensorField> vstFlds;
|
|
ReadFields(mesh, objects, vstFlds);
|
|
|
|
PtrList<volSymmTensorField> vsymtFlds;
|
|
ReadFields(mesh, objects, vsymtFlds);
|
|
|
|
PtrList<volTensorField> vtFlds;
|
|
ReadFields(mesh, objects, vtFlds);
|
|
|
|
// Read surface fields.
|
|
|
|
PtrList<surfaceScalarField> ssFlds;
|
|
ReadFields(mesh, objects, ssFlds);
|
|
|
|
PtrList<surfaceVectorField> svFlds;
|
|
ReadFields(mesh, objects, svFlds);
|
|
|
|
PtrList<surfaceSphericalTensorField> sstFlds;
|
|
ReadFields(mesh, objects, sstFlds);
|
|
|
|
PtrList<surfaceSymmTensorField> ssymtFlds;
|
|
ReadFields(mesh, objects, ssymtFlds);
|
|
|
|
PtrList<surfaceTensorField> stFlds;
|
|
ReadFields(mesh, objects, stFlds);
|
|
|
|
Info<< endl;
|
|
|
|
|
|
// Remove any demand-driven fields ('S', 'V' etc)
|
|
mesh.clearOut();
|
|
|
|
|
|
if (nCellRegions == 1)
|
|
{
|
|
Info<< "Only one region. Doing nothing." << endl;
|
|
}
|
|
else if (makeCellZones)
|
|
{
|
|
Info<< "Putting cells into cellZones instead of splitting mesh."
|
|
<< endl;
|
|
|
|
// Check if region overlaps with existing zone. If so keep.
|
|
|
|
for (label regionI = 0; regionI < nCellRegions; regionI++)
|
|
{
|
|
const labelList& zones = regionToZones[regionI];
|
|
|
|
if (zones.size() == 1 && zones[0] != -1)
|
|
{
|
|
// Existing zone
|
|
const label zoneI = zones[0];
|
|
Info<< " Region " << regionI << " : corresponds to existing"
|
|
<< " cellZone "
|
|
<< zoneI << ' ' << mesh.cellZones()[zoneI].name() << endl;
|
|
}
|
|
else
|
|
{
|
|
// Create new cellZone.
|
|
const labelList regionCells(findIndices(cellRegion, regionI));
|
|
|
|
const word zoneName = "region" + Foam::name(regionI);
|
|
|
|
label zoneI = mesh.cellZones().findZoneID(zoneName);
|
|
|
|
if (zoneI == -1)
|
|
{
|
|
zoneI = mesh.cellZones().size();
|
|
mesh.cellZones().setSize(zoneI+1);
|
|
mesh.cellZones().set
|
|
(
|
|
zoneI,
|
|
new cellZone
|
|
(
|
|
zoneName, //name
|
|
std::move(regionCells), //addressing
|
|
zoneI, //index
|
|
mesh.cellZones() //cellZoneMesh
|
|
)
|
|
);
|
|
}
|
|
else
|
|
{
|
|
mesh.cellZones()[zoneI].clearAddressing();
|
|
mesh.cellZones()[zoneI] = regionCells;
|
|
}
|
|
Info<< " Region " << regionI << " : created new cellZone "
|
|
<< zoneI << ' ' << mesh.cellZones()[zoneI].name() << endl;
|
|
}
|
|
}
|
|
mesh.cellZones().writeOpt(IOobject::AUTO_WRITE);
|
|
|
|
if (!overwrite)
|
|
{
|
|
++runTime;
|
|
mesh.setInstance(runTime.timeName());
|
|
}
|
|
else
|
|
{
|
|
mesh.setInstance(oldInstance);
|
|
}
|
|
|
|
Info<< "Writing cellZones as new mesh to time " << runTime.timeName()
|
|
<< nl << endl;
|
|
|
|
mesh.write();
|
|
|
|
|
|
// Write cellSets for convenience
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
Info<< "Writing cellSets corresponding to cellZones." << nl << endl;
|
|
|
|
for (const auto& cz : mesh.cellZones())
|
|
{
|
|
cellSet(mesh, cz.name(), cz).write();
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Add patches for interfaces
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
// Add all possible patches. Empty ones get filtered later on.
|
|
Info<< nl << "Adding patches" << nl << endl;
|
|
|
|
labelList interfacePatches
|
|
(
|
|
addRegionPatches
|
|
(
|
|
mesh,
|
|
regionNames,
|
|
interfaces,
|
|
interfaceNames
|
|
)
|
|
);
|
|
|
|
|
|
if (!overwrite)
|
|
{
|
|
++runTime;
|
|
}
|
|
|
|
|
|
// Create regions
|
|
// ~~~~~~~~~~~~~~
|
|
|
|
if (insidePoint)
|
|
{
|
|
const point insidePoint = args.get<point>("insidePoint");
|
|
|
|
label regionI = -1;
|
|
|
|
(void)mesh.tetBasePtIs();
|
|
|
|
label celli = mesh.findCell(insidePoint);
|
|
|
|
Info<< nl << "Found point " << insidePoint << " in cell " << celli
|
|
<< endl;
|
|
|
|
if (celli != -1)
|
|
{
|
|
regionI = cellRegion[celli];
|
|
}
|
|
|
|
reduce(regionI, maxOp<label>());
|
|
|
|
Info<< nl
|
|
<< "Subsetting region " << regionI
|
|
<< " containing point " << insidePoint << endl;
|
|
|
|
if (regionI == -1)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Point " << insidePoint
|
|
<< " is not inside the mesh." << nl
|
|
<< "Bounding box of the mesh:" << mesh.bounds()
|
|
<< exit(FatalError);
|
|
}
|
|
|
|
createAndWriteRegion
|
|
(
|
|
mesh,
|
|
cellRegion,
|
|
regionNames,
|
|
prefixRegion,
|
|
faceToInterface,
|
|
interfacePatches,
|
|
regionI,
|
|
(overwrite ? oldInstance : runTime.timeName())
|
|
);
|
|
}
|
|
else if (largestOnly)
|
|
{
|
|
label regionI = findMax(regionSizes);
|
|
|
|
Info<< nl
|
|
<< "Subsetting region " << regionI
|
|
<< " of size " << regionSizes[regionI]
|
|
<< " as named region " << regionNames[regionI] << endl;
|
|
|
|
createAndWriteRegion
|
|
(
|
|
mesh,
|
|
cellRegion,
|
|
regionNames,
|
|
prefixRegion,
|
|
faceToInterface,
|
|
interfacePatches,
|
|
regionI,
|
|
(overwrite ? oldInstance : runTime.timeName())
|
|
);
|
|
}
|
|
else
|
|
{
|
|
// Split all
|
|
for (label regionI = 0; regionI < nCellRegions; regionI++)
|
|
{
|
|
Info<< nl
|
|
<< "Region " << regionI << nl
|
|
<< "-------- " << endl;
|
|
|
|
createAndWriteRegion
|
|
(
|
|
mesh,
|
|
cellRegion,
|
|
regionNames,
|
|
prefixRegion,
|
|
faceToInterface,
|
|
interfacePatches,
|
|
regionI,
|
|
(overwrite ? oldInstance : runTime.timeName())
|
|
);
|
|
}
|
|
}
|
|
}
|
|
|
|
Info<< "End\n" << endl;
|
|
|
|
return 0;
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|