Rewrite + region addressing

This commit is contained in:
Hrvoje Jasak 2011-05-27 11:38:05 +01:00
parent 7c464b2a2e
commit 6a7971f2e4

View file

@ -31,7 +31,7 @@ Description
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).
- additional faces from faceSet (-blockedFaces faceSet).
- any face in between differing cellZones (-cellZones)
Output is:
@ -105,7 +105,8 @@ void addPatchFields(fvMesh& mesh, const word& patchFieldType)
for
(
typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
typename HashTable<const GeoField*>::const_iterator iter =
flds.begin();
iter != flds.end();
++iter
)
@ -145,7 +146,8 @@ void trimPatchFields(fvMesh& mesh, const label nPatches)
for
(
typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
typename HashTable<const GeoField*>::const_iterator iter =
flds.begin();
iter != flds.end();
++iter
)
@ -171,7 +173,8 @@ void reorderPatchFields(fvMesh& mesh, const labelList& oldToNew)
for
(
typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
typename HashTable<const GeoField*>::const_iterator iter =
flds.begin();
iter != flds.end();
++iter
)
@ -196,6 +199,7 @@ label addPatch(fvMesh& mesh, const polyPatch& patch)
const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
label patchI = polyPatches.findPatchID(patch.name());
if (patchI != -1)
{
if (polyPatches[patchI].type() == patch.type())
@ -239,6 +243,9 @@ label addPatch(fvMesh& mesh, const polyPatch& patch)
fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
Info<< "Inserting patch " << patch.name() << " to slot " << insertPatchI
<< " out of " << polyPatches.size() << endl;
// Add polyPatch at the end
polyPatches.setSize(sz + 1);
polyPatches.set
@ -351,47 +358,7 @@ label addPatch(fvMesh& mesh, const polyPatch& patch)
}
// Reorder and delete patches.
void reorderPatches
(
fvMesh& mesh,
const labelList& oldToNew,
const label nNewPatches
)
{
polyBoundaryMesh& polyPatches =
const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
// Shuffle into place
polyPatches.reorder(oldToNew);
fvPatches.reorder(oldToNew);
reorderPatchFields<volScalarField>(mesh, oldToNew);
reorderPatchFields<volVectorField>(mesh, oldToNew);
reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
reorderPatchFields<volTensorField>(mesh, oldToNew);
reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
// Remove last.
polyPatches.setSize(nNewPatches);
fvPatches.setSize(nNewPatches);
trimPatchFields<volScalarField>(mesh, nNewPatches);
trimPatchFields<volVectorField>(mesh, nNewPatches);
trimPatchFields<volSphericalTensorField>(mesh, nNewPatches);
trimPatchFields<volSymmTensorField>(mesh, nNewPatches);
trimPatchFields<volTensorField>(mesh, nNewPatches);
trimPatchFields<surfaceScalarField>(mesh, nNewPatches);
trimPatchFields<surfaceVectorField>(mesh, nNewPatches);
trimPatchFields<surfaceSphericalTensorField>(mesh, nNewPatches);
trimPatchFields<surfaceSymmTensorField>(mesh, nNewPatches);
trimPatchFields<surfaceTensorField>(mesh, nNewPatches);
}
// Reorder and delete patches. Deleted. HJ, 22/May/2011
template<class GeoField>
@ -400,10 +367,12 @@ void subsetVolFields
const fvMesh& mesh,
const fvMesh& subMesh,
const labelList& cellMap,
const labelList& faceMap
const labelList& faceMap,
const labelList& patchMap
)
{
const labelList patchMap(identity(mesh.boundaryMesh().size()));
// Real patch map passed by argument
// HJ, 22/May/2011
HashTable<const GeoField*> fields
(
@ -459,15 +428,18 @@ void subsetSurfaceFields
(
const fvMesh& mesh,
const fvMesh& subMesh,
const labelList& faceMap
const labelList& faceMap,
const labelList& patchMap
)
{
const labelList patchMap(identity(mesh.boundaryMesh().size()));
// Real patch map passed by argument
// HJ, 22/May/2011
HashTable<const GeoField*> fields
(
mesh.objectRegistry::lookupClass<GeoField>()
);
forAllConstIter(typename HashTable<const GeoField*>, fields, iter)
{
const GeoField& fld = *iter();
@ -511,6 +483,7 @@ void subsetSurfaceFields
}
}
// Select all cells not in the region
labelList getNonRegionCells(const labelList& cellRegion, const label regionI)
{
@ -678,69 +651,19 @@ void getInterfaceSizes
}
// Create mesh for region.
// Create mesh for region. Mesh pointer not passed in: instead, mesh is
// written out (as polyMesh) and read in (as fvMesh)
// to allow for field mapping.
// HJ, 5/Apr/2011
autoPtr<mapPolyMesh> createRegionMesh
(
const labelList& cellRegion,
const EdgeMap<label>& interfaceToPatch,
const fvMesh& mesh,
const label regionI,
const word& regionName,
autoPtr<fvMesh>& newMesh
const word& regionName
)
{
// Create dummy system/fv*
{
IOobject io
(
"fvSchemes",
mesh.time().system(),
regionName,
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
);
Info<< "Testing:" << io.objectPath() << endl;
if (!io.headerOk())
// if (!exists(io.objectPath()))
{
Info<< "Writing dummy " << regionName/io.name() << endl;
dictionary dummyDict;
dictionary divDict;
dummyDict.add("divSchemes", divDict);
dictionary gradDict;
dummyDict.add("gradSchemes", gradDict);
dictionary laplDict;
dummyDict.add("laplacianSchemes", laplDict);
IOdictionary(io, dummyDict).regIOobject::write();
}
}
{
IOobject io
(
"fvSolution",
mesh.time().system(),
regionName,
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
);
if (!io.headerOk())
//if (!exists(io.objectPath()))
{
Info<< "Writing dummy " << regionName/io.name() << endl;
dictionary dummyDict;
IOdictionary(io, dummyDict).regIOobject::write();
}
}
// Neighbour cellRegion.
labelList coupledRegion(mesh.nFaces() - mesh.nInternalFaces());
@ -751,7 +674,6 @@ autoPtr<mapPolyMesh> createRegionMesh
}
syncTools::swapBoundaryFaceList(mesh, coupledRegion, false);
// Topology change container. Start off from existing mesh.
directTopoChange meshMod(mesh);
@ -829,6 +751,8 @@ autoPtr<mapPolyMesh> createRegionMesh
meshMod
);
autoPtr<polyMesh> newMesh;
autoPtr<mapPolyMesh> map = meshMod.makeMesh
(
newMesh,
@ -843,109 +767,8 @@ autoPtr<mapPolyMesh> createRegionMesh
mesh
);
return map;
}
void createAndWriteRegion
(
const fvMesh& mesh,
const labelList& cellRegion,
const wordList& regionNames,
const EdgeMap<label>& interfaceToPatch,
const label regionI,
const word& newMeshInstance
)
{
Info<< "Creating mesh for region " << regionI
<< ' ' << regionNames[regionI] << endl;
autoPtr<fvMesh> newMesh;
autoPtr<mapPolyMesh> map = createRegionMesh
(
cellRegion,
interfaceToPatch,
mesh,
regionI,
regionNames[regionI],
newMesh
);
Info<< "Mapping fields" << endl;
// Map existing fields
newMesh().updateMesh(map());
// Add subsetted fields
subsetVolFields<volScalarField>
(
mesh,
newMesh(),
map().cellMap(),
map().faceMap()
);
subsetVolFields<volVectorField>
(
mesh,
newMesh(),
map().cellMap(),
map().faceMap()
);
subsetVolFields<volSphericalTensorField>
(
mesh,
newMesh(),
map().cellMap(),
map().faceMap()
);
subsetVolFields<volSymmTensorField>
(
mesh,
newMesh(),
map().cellMap(),
map().faceMap()
);
subsetVolFields<volTensorField>
(
mesh,
newMesh(),
map().cellMap(),
map().faceMap()
);
subsetSurfaceFields<surfaceScalarField>
(
mesh,
newMesh(),
map().faceMap()
);
subsetSurfaceFields<surfaceVectorField>
(
mesh,
newMesh(),
map().faceMap()
);
subsetSurfaceFields<surfaceSphericalTensorField>
(
mesh,
newMesh(),
map().faceMap()
);
subsetSurfaceFields<surfaceSymmTensorField>
(
mesh,
newMesh(),
map().faceMap()
);
subsetSurfaceFields<surfaceTensorField>
(
mesh,
newMesh(),
map().faceMap()
);
const polyBoundaryMesh& newPatches = newMesh().boundaryMesh();
polyBoundaryMesh& newPatches =
const_cast<polyBoundaryMesh&>(newMesh().boundaryMesh());
newPatches.checkParallelSync(true);
// Delete empty patches
@ -993,18 +816,22 @@ void createAndWriteRegion
}
}
reorderPatches(newMesh(), oldToNew, nNewPatches);
polyBoundaryMesh& polyPatches =
const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
// Reorder polyPatches and trim empty patches from the end
// of the list
newPatches.reorder(oldToNew);
newPatches.setSize(nNewPatches);
// Write the mesh
Info<< "Writing new mesh" << endl;
newMesh().setInstance(newMeshInstance);
newMesh().write();
// Write addressing files like decomposePar
Info<< "Writing addressing to base mesh" << endl;
labelIOList pointProcAddressing
labelIOList pointRegionAddressing
(
IOobject
(
@ -1018,12 +845,13 @@ void createAndWriteRegion
),
map().pointMap()
);
Info<< "Writing map " << pointProcAddressing.name()
Info<< "Writing map " << pointRegionAddressing.name()
<< " from region" << regionI
<< " points back to base mesh." << endl;
pointProcAddressing.write();
pointRegionAddressing.write();
labelIOList faceProcAddressing
labelIOList faceRegionAddressing
(
IOobject
(
@ -1037,7 +865,8 @@ void createAndWriteRegion
),
newMesh().nFaces()
);
forAll(faceProcAddressing, faceI)
forAll (faceRegionAddressing, faceI)
{
// face + turning index. (see decomposePar)
// Is the face pointing in the same direction?
@ -1049,19 +878,20 @@ void createAndWriteRegion
== mesh.faceOwner()[oldFaceI]
)
{
faceProcAddressing[faceI] = oldFaceI+1;
faceRegionAddressing[faceI] = oldFaceI + 1;
}
else
{
faceProcAddressing[faceI] = -(oldFaceI+1);
faceRegionAddressing[faceI] = -(oldFaceI + 1);
}
}
Info<< "Writing map " << faceProcAddressing.name()
Info<< "Writing map " << faceRegionAddressing.name()
<< " from region" << regionI
<< " faces back to base mesh." << endl;
faceProcAddressing.write();
faceRegionAddressing.write();
labelIOList cellProcAddressing
labelIOList cellRegionAddressing
(
IOobject
(
@ -1075,10 +905,197 @@ void createAndWriteRegion
),
map().cellMap()
);
Info<< "Writing map " <<cellProcAddressing.name()
Info<< "Writing map " << cellRegionAddressing.name()
<< " from region" << regionI
<< " cells back to base mesh." << endl;
cellProcAddressing.write();
cellRegionAddressing.write();
// Calculate and write boundary addressing
{
const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
const polyBoundaryMesh& newPatches = newMesh().boundaryMesh();
labelIOList boundaryRegionAddressing
(
IOobject
(
"boundaryRegionAddressing",
newMesh().facesInstance(),
newMesh().meshSubDir,
newMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
labelList
(
newMesh().boundaryMesh().size(),
-1
)
);
Info<<"Number of new patches: " << nNewPatches <<endl;
forAll (boundaryRegionAddressing, patchID)
{
const word& curPatchName = newPatches[patchID].name();
forAll (oldPatches, oldPatchI)
{
if (oldPatches[oldPatchI].name() == curPatchName)
{
boundaryRegionAddressing[patchID] = oldPatchI;
}
}
}
Info<< "Writing map " << boundaryRegionAddressing.name()
<< " from region" << regionI
<< " faces back to base mesh." << endl;
boundaryRegionAddressing.write();
}
return map;
}
void createAndWriteRegion
(
const fvMesh& mesh,
const labelList& cellRegion,
const wordList& regionNames,
const EdgeMap<label>& interfaceToPatch,
const label regionI,
const word& newMeshInstance
)
{
Info<< "Creating mesh for region " << regionI
<< ' ' << regionNames[regionI] << endl;
// Create region mesh will write the mesh
// HJ, 19/May/2011
autoPtr<mapPolyMesh> map = createRegionMesh
(
cellRegion,
interfaceToPatch,
mesh,
regionI,
regionNames[regionI]
);
// Read in mesh as fvMesh for mapping. HJ, 19/May/2011
fvMesh newMesh
(
IOobject
(
regionNames[regionI],
newMeshInstance,
mesh.time(),
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
// Read in patch map. HJ, 22/May/2011
labelIOList boundaryRegionAddressing
(
IOobject
(
"boundaryRegionAddressing",
newMesh.facesInstance(),
newMesh.meshSubDir,
newMesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
Info<< "Mapping fields" << endl;
// Map existing fields: not needed
// Add subset fields
// Note: real patch map provided in place of identity map hack
// HJ, 22/May/2011
subsetVolFields<volScalarField>
(
mesh,
newMesh,
map().cellMap(),
map().faceMap(),
boundaryRegionAddressing
);
subsetVolFields<volVectorField>
(
mesh,
newMesh,
map().cellMap(),
map().faceMap(),
boundaryRegionAddressing
);
subsetVolFields<volSphericalTensorField>
(
mesh,
newMesh,
map().cellMap(),
map().faceMap(),
boundaryRegionAddressing
);
subsetVolFields<volSymmTensorField>
(
mesh,
newMesh,
map().cellMap(),
map().faceMap(),
boundaryRegionAddressing
);
subsetVolFields<volTensorField>
(
mesh,
newMesh,
map().cellMap(),
map().faceMap(),
boundaryRegionAddressing
);
subsetSurfaceFields<surfaceScalarField>
(
mesh,
newMesh,
map().faceMap(),
boundaryRegionAddressing
);
subsetSurfaceFields<surfaceVectorField>
(
mesh,
newMesh,
map().faceMap(),
boundaryRegionAddressing
);
subsetSurfaceFields<surfaceSphericalTensorField>
(
mesh,
newMesh,
map().faceMap(),
boundaryRegionAddressing
);
subsetSurfaceFields<surfaceSymmTensorField>
(
mesh,
newMesh,
map().faceMap(),
boundaryRegionAddressing
);
subsetSurfaceFields<surfaceTensorField>
(
mesh,
newMesh,
map().faceMap(),
boundaryRegionAddressing
);
// Moved map writing into createRegionMesh. HJ, 20/May/2011
}
@ -1149,6 +1166,7 @@ EdgeMap<label> addRegionPatches
interfaceToPatch.insert(e, patchI);
}
}
return interfaceToPatch;
}
@ -1207,72 +1225,6 @@ label findCorrespondingRegion
}
//// Checks if cellZone has corresponding cellRegion.
//label findCorrespondingRegion
//(
// const cellZoneMesh& cellZones,
// const labelList& existingZoneID, // per cell the (unique) zoneID
// const labelList& cellRegion,
// const label nCellRegions,
// const label zoneI
//)
//{
// // Region corresponding to zone. Start off with special value: no
// // corresponding region.
// label regionI = labelMax;
//
// const cellZone& cz = cellZones[zoneI];
//
// if (cz.empty())
// {
// // My local portion is empty. Maps to any empty cellZone. Mark with
// // special value which can get overwritten by other processors.
// regionI = -1;
// }
// else
// {
// regionI = cellRegion[cz[0]];
//
// forAll(cz, i)
// {
// if (cellRegion[cz[i]] != regionI)
// {
// regionI = labelMax;
// break;
// }
// }
// }
//
// // Determine same zone over all processors.
// reduce(regionI, maxOp<label>());
//
//
// // 2. All of region present?
//
// if (regionI == labelMax)
// {
// regionI = -1;
// }
// else if (regionI != -1)
// {
// 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;
//}
// Get zone per cell
// - non-unique zoning
// - coupled zones
@ -1463,6 +1415,7 @@ void writeCellToRegion(const fvMesh& mesh, const labelList& cellRegion)
Info<< "Writing region per cell file (for manual decomposition) to "
<< cellToRegion.objectPath() << nl << endl;
}
// Write for postprocessing
{
volScalarField cellToRegion
@ -1549,8 +1502,6 @@ int main(int argc, char *argv[])
<< exit(FatalError);
}
if (insidePoint && largestOnly)
{
FatalErrorIn(args.executable())
@ -1565,8 +1516,10 @@ int main(int argc, char *argv[])
// Existing zoneID
labelList zoneID(mesh.nCells(), -1);
// Neighbour zoneID.
// Neighbour zoneID
labelList neiZoneID(mesh.nFaces() - mesh.nInternalFaces());
getZoneID(mesh, cellZones, zoneID, neiZoneID);
@ -1576,10 +1529,13 @@ int main(int argc, char *argv[])
// 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;
@ -1591,6 +1547,7 @@ int main(int argc, char *argv[])
<< " cells to be in one and only one cellZone." << nl << endl;
label unzonedCellI = findIndex(zoneID, -1);
if (unzonedCellI != -1)
{
FatalErrorIn(args.executable())
@ -1601,11 +1558,13 @@ int main(int argc, char *argv[])
<< " 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;
@ -1640,6 +1599,7 @@ int main(int argc, char *argv[])
getZoneID(mesh, newCellZones, newZoneID, newNeiZoneID);
label unzonedCellI = findIndex(newZoneID, -1);
if (unzonedCellI != -1)
{
FatalErrorIn(args.executable())
@ -1650,11 +1610,13 @@ int main(int argc, char *argv[])
<< " 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;
@ -1676,7 +1638,8 @@ int main(int argc, char *argv[])
faceSet blockedFaceSet(mesh, blockedFacesName);
Info<< "Read "
<< returnReduce(blockedFaceSet.size(), sumOp<label>())
<< " blocked faces from set " << blockedFacesName << nl << endl;
<< " blocked faces from set " << blockedFacesName
<< nl << endl;
blockedFace.setSize(mesh.nFaces(), false);
@ -1716,6 +1679,7 @@ int main(int argc, char *argv[])
// Do a topological walk to determine regions
regionSplit regions(mesh, blockedFace);
nCellRegions = regions.nRegions();
cellRegion.transfer(regions);
@ -1740,7 +1704,6 @@ int main(int argc, char *argv[])
writeCellToRegion(mesh, cellRegion);
// Sizes per region
// ~~~~~~~~~~~~~~~~