Merge commit 'OpenCFD/master' into olesenm

This commit is contained in:
Mark Olesen 2009-09-23 09:14:28 +02:00
commit ad5d10b817
16 changed files with 543 additions and 189 deletions

View File

@ -64,7 +64,7 @@ class calcEntry
public:
//- Runtime type information
TypeName("calc");
ClassName("calc");
// Member Functions

View File

@ -32,7 +32,8 @@ Description
patch_YYY_XXX.obj : all face centres of patch YYY
Optional: patch faces (as polygons) : patchFaces_YYY_XXX.obj
Optional: - patch faces (as polygons) : patchFaces_YYY_XXX.obj
- non-manifold edges : patchEdges_YYY_XXX.obj
\*---------------------------------------------------------------------------*/
@ -51,7 +52,7 @@ using namespace Foam;
void writeOBJ(const point& pt, Ostream& os)
{
os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
}
// All edges of mesh
@ -75,8 +76,7 @@ void writePoints(const polyMesh& mesh, const fileName& timeName)
{
const edge& e = mesh.edges()[edgeI];
pointStream << "l " << e.start() + 1 << ' ' << e.end() + 1
<< endl;
pointStream << "l " << e.start() + 1 << ' ' << e.end() + 1 << nl;
}
}
@ -277,7 +277,47 @@ void writePatchFaces
{
patchFaceStream << ' ' << f[fp]+1;
}
patchFaceStream << endl;
patchFaceStream << nl;
}
}
}
void writePatchBoundaryEdges
(
const polyMesh& mesh,
const fileName& timeName
)
{
const polyBoundaryMesh& patches = mesh.boundaryMesh();
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
fileName edgeFile
(
mesh.time().path()
/ "patchEdges_" + pp.name() + '_' + timeName + ".obj"
);
Info << "Writing patch edges to " << edgeFile << endl;
OFstream patchEdgeStream(edgeFile);
forAll(pp.localPoints(), pointI)
{
writeOBJ(pp.localPoints()[pointI], patchEdgeStream);
}
for (label edgeI = pp.nInternalEdges(); edgeI < pp.nEdges(); edgeI++)
{
if (pp.edgeFaces()[edgeI].size() == 1)
{
const edge& e = pp.edges()[edgeI];
patchEdgeStream<< "l " << e[0]+1 << ' ' << e[1]+1 << nl;
}
}
}
}
@ -339,6 +379,7 @@ int main(int argc, char *argv[])
{
timeSelector::addOptions();
argList::validOptions.insert("patchFaces", "");
argList::validOptions.insert("patchEdges", "");
argList::validOptions.insert("cell", "cellI");
argList::validOptions.insert("face", "faceI");
argList::validOptions.insert("point", "pointI");
@ -351,6 +392,7 @@ int main(int argc, char *argv[])
runTime.functionObjects().off();
bool patchFaces = args.optionFound("patchFaces");
bool patchEdges = args.optionFound("patchEdges");
bool doCell = args.optionFound("cell");
bool doPoint = args.optionFound("point");
bool doFace = args.optionFound("face");
@ -381,19 +423,23 @@ int main(int argc, char *argv[])
{
writePatchFaces(mesh, runTime.timeName());
}
else if (doCell)
if (patchEdges)
{
writePatchBoundaryEdges(mesh, runTime.timeName());
}
if (doCell)
{
label cellI = args.optionRead<label>("cell");
writePoints(mesh, cellI, runTime.timeName());
}
else if (doPoint)
if (doPoint)
{
label pointI = args.optionRead<label>("point");
writePointCells(mesh, pointI, runTime.timeName());
}
else if (doFace)
if (doFace)
{
label faceI = args.optionRead<label>("face");
@ -415,7 +461,7 @@ int main(int argc, char *argv[])
meshTools::writeOBJ(str, faceList(1, f), mesh.points());
}
else if (doCellSet)
if (doCellSet)
{
word setName(args.option("cellSet"));
@ -427,7 +473,7 @@ int main(int argc, char *argv[])
writePoints(mesh, cells.toc(), runTime.timeName());
}
else if (doFaceSet)
if (doFaceSet)
{
word setName(args.option("faceSet"));
@ -458,7 +504,16 @@ int main(int argc, char *argv[])
faces.toc()
);
}
else
else if
(
!patchFaces
&& !patchEdges
&& !doCell
&& !doPoint
&& !doFace
&& !doCellSet
&& !doFaceSet
)
{
// points & edges
writePoints(mesh, runTime.timeName());

View File

@ -71,7 +71,8 @@ linearNormalCoeffs
linearDirectionCoeffs
{
direction (0 0 1);
direction (0 1 0);
thickness 0.05;
}
linearRadialCoeffs
@ -89,7 +90,7 @@ sigmaRadialCoeffs
// Do front and back need to be merged? Usually only makes sense for 360
// degree wedges.
mergeFaces true;
mergeFaces false; //true;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -49,7 +49,7 @@ int main(int argc, char *argv[])
(
IOobject
(
mergePolyMesh::defaultRegion,
masterRegion,
runTimeMaster.timeName(),
runTimeMaster
)
@ -63,7 +63,7 @@ int main(int argc, char *argv[])
(
IOobject
(
mergePolyMesh::defaultRegion,
addRegion,
runTimeToAdd.timeName(),
runTimeToAdd
)

View File

@ -2,9 +2,11 @@
argList::validArgs.append("master root");
argList::validArgs.append("master case");
argList::validOptions.insert("masterRegion", "name");
argList::validArgs.append("root to add");
argList::validArgs.append("case to add");
argList::validOptions.insert("addRegion", "name");
argList args(argc, argv);
@ -15,9 +17,15 @@
fileName rootDirMaster(args.additionalArgs()[0]);
fileName caseDirMaster(args.additionalArgs()[1]);
word masterRegion = polyMesh::defaultRegion;
args.optionReadIfPresent("masterRegion", masterRegion);
fileName rootDirToAdd(args.additionalArgs()[2]);
fileName caseDirToAdd(args.additionalArgs()[3]);
word addRegion = polyMesh::defaultRegion;
args.optionReadIfPresent("addRegion", addRegion);
Info<< "Master: " << rootDirMaster << " " << caseDirMaster << nl
<< "mesh to add: " << rootDirToAdd << " " << caseDirToAdd << endl;
Info<< "Master: " << rootDirMaster << " " << caseDirMaster
<< " region " << masterRegion << nl
<< "mesh to add: " << rootDirToAdd << " " << caseDirToAdd
<< " region " << addRegion << endl;

View File

@ -43,6 +43,9 @@ Description
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) region.
\*---------------------------------------------------------------------------*/
@ -609,7 +612,7 @@ void getInterfaceSizes
// Create mesh for region.
autoPtr<mapPolyMesh> createRegionMesh
(
const regionSplit& cellRegion,
const labelList& cellRegion,
const EdgeMap<label>& interfaceToPatch,
const fvMesh& mesh,
const label regionI,
@ -766,7 +769,7 @@ autoPtr<mapPolyMesh> createRegionMesh
void createAndWriteRegion
(
const fvMesh& mesh,
const regionSplit& cellRegion,
const labelList& cellRegion,
const wordList& regionNames,
const EdgeMap<label>& interfaceToPatch,
const label regionI,
@ -1005,7 +1008,8 @@ void createAndWriteRegion
EdgeMap<label> addRegionPatches
(
fvMesh& mesh,
const regionSplit& cellRegion,
const labelList& cellRegion,
const label nCellRegions,
const edgeList& interfaces,
const EdgeMap<label>& interfaceSizes,
const wordList& regionNames
@ -1016,7 +1020,7 @@ EdgeMap<label> addRegionPatches
Info<< nl << "Adding patches" << nl << endl;
EdgeMap<label> interfaceToPatch(cellRegion.nRegions());
EdgeMap<label> interfaceToPatch(nCellRegions);
forAll(interfaces, interI)
{
@ -1108,13 +1112,14 @@ EdgeMap<label> addRegionPatches
label findCorrespondingRegion
(
const labelList& existingZoneID, // per cell the (unique) zoneID
const regionSplit& cellRegion,
const labelList& cellRegion,
const label nCellRegions,
const label zoneI,
const label minOverlapSize
)
{
// Per region the number of cells in zoneI
labelList cellsInZone(cellRegion.nRegions(), 0);
labelList cellsInZone(nCellRegions, 0);
forAll(cellRegion, cellI)
{
@ -1162,7 +1167,8 @@ label findCorrespondingRegion
//(
// const cellZoneMesh& cellZones,
// const labelList& existingZoneID, // per cell the (unique) zoneID
// const regionSplit& cellRegion,
// const labelList& cellRegion,
// const label nCellRegions,
// const label zoneI
//)
//{
@ -1227,6 +1233,7 @@ label findCorrespondingRegion
int main(int argc, char *argv[])
{
argList::validOptions.insert("cellZones", "");
argList::validOptions.insert("cellZonesOnly", "");
argList::validOptions.insert("blockedFaces", "faceSet");
argList::validOptions.insert("makeCellZones", "");
argList::validOptions.insert("largestOnly", "");
@ -1249,13 +1256,14 @@ int main(int argc, char *argv[])
<< blockedFacesName << nl << endl;
}
bool makeCellZones = args.optionFound("makeCellZones");
bool largestOnly = args.optionFound("largestOnly");
bool insidePoint = args.optionFound("insidePoint");
bool useCellZones = args.optionFound("cellZones");
bool overwrite = args.optionFound("overwrite");
bool detectOnly = args.optionFound("detectOnly");
bool sloppyCellZones = args.optionFound("sloppyCellZones");
bool makeCellZones = args.optionFound("makeCellZones");
bool largestOnly = args.optionFound("largestOnly");
bool insidePoint = args.optionFound("insidePoint");
bool useCellZones = args.optionFound("cellZones");
bool useCellZonesOnly = args.optionFound("cellZonesOnly");
bool overwrite = args.optionFound("overwrite");
bool detectOnly = args.optionFound("detectOnly");
bool sloppyCellZones = args.optionFound("sloppyCellZones");
if (insidePoint && largestOnly)
{
@ -1370,13 +1378,37 @@ int main(int argc, char *argv[])
}
// Do the topological walk to determine regions
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Determine per cell the region it belongs to
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// regionSplit is the labelList with the region per cell.
regionSplit cellRegion(mesh, blockedFace);
// 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:" << cellRegion.nRegions() << nl << endl;
Info<< endl << "Number of regions:" << nCellRegions << nl << endl;
// Write to manual decomposition option
@ -1429,7 +1461,7 @@ int main(int argc, char *argv[])
// Sizes per region
// ~~~~~~~~~~~~~~~~
labelList regionSizes(cellRegion.nRegions(), 0);
labelList regionSizes(nCellRegions, 0);
forAll(cellRegion, cellI)
{
@ -1489,9 +1521,9 @@ int main(int argc, char *argv[])
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Region per zone
labelList regionToZone(cellRegion.nRegions(), -1);
labelList regionToZone(nCellRegions, -1);
// Name of region
wordList regionNames(cellRegion.nRegions());
wordList regionNames(nCellRegions);
// Zone to region
labelList zoneToRegion(cellZones.size(), -1);
@ -1506,6 +1538,7 @@ int main(int argc, char *argv[])
(
zoneID,
cellRegion,
nCellRegions,
zoneI,
label(0.5*zoneSizes[zoneI]) // minimum overlap
);
@ -1532,6 +1565,7 @@ int main(int argc, char *argv[])
(
zoneID,
cellRegion,
nCellRegions,
zoneI,
1 // minimum overlap
);
@ -1660,7 +1694,7 @@ int main(int argc, char *argv[])
mesh.clearOut();
if (cellRegion.nRegions() == 1)
if (nCellRegions == 1)
{
Info<< "Only one region. Doing nothing." << endl;
}
@ -1671,7 +1705,7 @@ int main(int argc, char *argv[])
// Check if region overlaps with existing zone. If so keep.
for (label regionI = 0; regionI < cellRegion.nRegions(); regionI++)
for (label regionI = 0; regionI < nCellRegions; regionI++)
{
label zoneI = regionToZone[regionI];
@ -1759,6 +1793,7 @@ int main(int argc, char *argv[])
(
mesh,
cellRegion,
nCellRegions,
interfaces,
interfaceSizes,
regionNames
@ -1837,7 +1872,7 @@ int main(int argc, char *argv[])
else
{
// Split all
for (label regionI = 0; regionI < cellRegion.nRegions(); regionI++)
for (label regionI = 0; regionI < nCellRegions; regionI++)
{
Info<< nl
<< "Region " << regionI << nl

View File

@ -38,13 +38,14 @@ Description
int main(int argc, char *argv[])
{
# include "addRegionOption.H"
timeSelector::addOptions();
argList::validArgs.append("fieldName");
argList::validArgs.append("patchName");
# include "setRootCase.H"
# include "createTime.H"
instantList timeDirs = timeSelector::select0(runTime, args);
# include "createMesh.H"
# include "createNamedMesh.H"
word fieldName(args.additionalArgs()[0]);
word patchName(args.additionalArgs()[1]);

View File

@ -81,7 +81,12 @@ bool Foam::primitiveEntry::expandVariable
word varName = w(1, w.size()-1);
// lookup the variable name in the given dictionary....
const entry* ePtr = dict.lookupEntryPtr(varName, true, true);
// Note: allow wildcards to match? For now disabled since following
// would expand internalField to wildcard match and not expected
// internalField:
// internalField XXX;
// boundaryField { ".*" {YYY;} movingWall {value $internalField;}
const entry* ePtr = dict.lookupEntryPtr(varName, true, false);
// ...if defined insert its tokens into this
if (ePtr != NULL)

View File

@ -2509,43 +2509,55 @@ void Foam::autoLayerDriver::addLayers
(
const layerParameters& layerParams,
const dictionary& motionDict,
const labelList& patchIDs,
const label nAllowableErrors,
motionSmoother& meshMover,
decompositionMethod& decomposer,
fvMeshDistribute& distributor
)
{
fvMesh& mesh = meshRefiner_.mesh();
const indirectPrimitivePatch& pp = meshMover.patch();
const labelList& meshPoints = pp.meshPoints();
// Precalculate mesh edge labels for patch edges
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
labelList meshEdges(pp.nEdges());
forAll(pp.edges(), edgeI)
{
const edge& ppEdge = pp.edges()[edgeI];
label v0 = meshPoints[ppEdge[0]];
label v1 = meshPoints[ppEdge[1]];
meshEdges[edgeI] = meshTools::findEdge
autoPtr<indirectPrimitivePatch> pp
(
meshRefinement::makePatch
(
mesh.edges(),
mesh.pointEdges()[v0],
v0,
v1
);
}
mesh,
patchIDs
)
);
// Construct iterative mesh mover.
Info<< "Constructing mesh displacer ..." << endl;
autoPtr<motionSmoother> meshMover
(
new motionSmoother
(
mesh,
pp(),
patchIDs,
meshRefinement::makeDisplacementField
(
pointMesh::New(mesh),
patchIDs
),
motionDict
)
);
// Point-wise extrusion data
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Displacement for all pp.localPoints.
vectorField patchDisp(pp.nPoints(), vector::one);
vectorField patchDisp(pp().nPoints(), vector::one);
// Number of layers for all pp.localPoints. Note: only valid if
// extrudeStatus = EXTRUDE.
labelList patchNLayers(pp.nPoints(), 0);
labelList patchNLayers(pp().nPoints(), 0);
// Whether to add edge for all pp.localPoints.
List<extrudeMode> extrudeStatus(pp.nPoints(), EXTRUDE);
List<extrudeMode> extrudeStatus(pp().nPoints(), EXTRUDE);
// Get number of layer per point from number of layers per patch
@ -2554,7 +2566,7 @@ void Foam::autoLayerDriver::addLayers
setNumLayers
(
layerParams.numLayers(), // per patch the num layers
meshMover.adaptPatchIDs(), // patches that are being moved
meshMover().adaptPatchIDs(),// patches that are being moved
pp, // indirectpatch for all faces moving
patchDisp,
@ -2562,6 +2574,9 @@ void Foam::autoLayerDriver::addLayers
extrudeStatus
);
// Precalculate mesh edge labels for patch edges
labelList meshEdges(pp().meshEdges(mesh.edges(), mesh.pointEdges()));
// Disable extrusion on non-manifold points
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -2582,7 +2597,7 @@ void Foam::autoLayerDriver::addLayers
(
pp,
meshEdges,
layerParams.featureAngle()*constant::math::pi/180.0,
layerParams.featureAngle()*constant::math::pi/180.0,
patchDisp,
patchNLayers,
@ -2633,14 +2648,15 @@ void Foam::autoLayerDriver::addLayers
);
}
// Determine (wanted) point-wise layer thickness and expansion ratio
scalarField thickness(pp.nPoints());
scalarField minThickness(pp.nPoints());
scalarField expansionRatio(pp.nPoints());
scalarField thickness(pp().nPoints());
scalarField minThickness(pp().nPoints());
scalarField expansionRatio(pp().nPoints());
calculateLayerThickness
(
pp,
meshMover.adaptPatchIDs(),
meshMover().adaptPatchIDs(),
layerParams.expansionRatio(),
layerParams.relativeSizes(), // thickness relative to cellsize?
@ -2663,11 +2679,11 @@ void Foam::autoLayerDriver::addLayers
// Find maximum length of a patch name, for a nicer output
label maxPatchNameLen = 0;
forAll(meshMover.adaptPatchIDs(), i)
forAll(meshMover().adaptPatchIDs(), i)
{
label patchI = meshMover.adaptPatchIDs()[i];
label patchI = meshMover().adaptPatchIDs()[i];
word patchName = patches[patchI].name();
maxPatchNameLen = max(maxPatchNameLen,label(patchName.size()));
maxPatchNameLen = max(maxPatchNameLen, label(patchName.size()));
}
Info<< nl
@ -2678,9 +2694,9 @@ void Foam::autoLayerDriver::addLayers
<< setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
<< setw(0) << " ----- ------ --------- -------" << endl;
forAll(meshMover.adaptPatchIDs(), i)
forAll(meshMover().adaptPatchIDs(), i)
{
label patchI = meshMover.adaptPatchIDs()[i];
label patchI = meshMover().adaptPatchIDs()[i];
const labelList& meshPoints = patches[patchI].meshPoints();
@ -2691,7 +2707,7 @@ void Foam::autoLayerDriver::addLayers
forAll(meshPoints, patchPointI)
{
label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
label ppPointI = pp().meshPointMap()[meshPoints[patchPointI]];
//maxThickness = max(maxThickness, thickness[ppPointI]);
//minThickness = min(minThickness, thickness[ppPointI]);
@ -2761,7 +2777,7 @@ void Foam::autoLayerDriver::addLayers
IOobject::NO_WRITE,
false
),
meshMover.pMesh(),
meshMover().pMesh(),
dimensionedScalar("pointMedialDist", dimless, 0.0)
);
@ -2776,7 +2792,7 @@ void Foam::autoLayerDriver::addLayers
IOobject::NO_WRITE,
false
),
meshMover.pMesh(),
meshMover().pMesh(),
dimensionedVector("dispVec", dimless, vector::zero)
);
@ -2791,7 +2807,7 @@ void Foam::autoLayerDriver::addLayers
IOobject::NO_WRITE,
false
),
meshMover.pMesh(),
meshMover().pMesh(),
dimensionedScalar("medialRatio", dimless, 0.0)
);
@ -2871,7 +2887,7 @@ void Foam::autoLayerDriver::addLayers
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
{
pointField oldPatchPos(pp.localPoints());
pointField oldPatchPos(pp().localPoints());
//// Laplacian displacement shrinking.
//shrinkMeshDistance
@ -2886,7 +2902,7 @@ void Foam::autoLayerDriver::addLayers
// Medial axis based shrinking
shrinkMeshMedialDistance
(
meshMover,
meshMover(),
meshQualityDict,
layerParams.nSmoothThickness(),
@ -2908,7 +2924,7 @@ void Foam::autoLayerDriver::addLayers
);
// Update patchDisp (since not all might have been honoured)
patchDisp = oldPatchPos - pp.localPoints();
patchDisp = oldPatchPos - pp().localPoints();
}
// Truncate displacements that are too small (this will do internal
@ -2916,7 +2932,7 @@ void Foam::autoLayerDriver::addLayers
faceSet dummySet(mesh, "wrongPatchFaces", 0);
truncateDisplacement
(
meshMover,
meshMover(),
minThickness,
dummySet,
patchDisp,
@ -2931,7 +2947,7 @@ void Foam::autoLayerDriver::addLayers
dumpDisplacement
(
mesh.time().path()/"layer",
pp,
pp(),
patchDisp,
extrudeStatus
);
@ -2955,11 +2971,11 @@ void Foam::autoLayerDriver::addLayers
// Determine per point/per face number of layers to extrude. Also
// handles the slow termination of layers when going switching layers
labelList nPatchPointLayers(pp.nPoints(),-1);
labelList nPatchFaceLayers(pp.localFaces().size(),-1);
labelList nPatchPointLayers(pp().nPoints(),-1);
labelList nPatchFaceLayers(pp().localFaces().size(),-1);
setupLayerInfoTruncation
(
meshMover,
meshMover(),
patchNLayers,
extrudeStatus,
layerParams.nBufferCellsNoExtrude(),
@ -2998,7 +3014,7 @@ void Foam::autoLayerDriver::addLayers
addLayer.setRefinement
(
invExpansionRatio,
pp,
pp(),
labelList(0), // exposed patchIDs, not used for adding layers
nPatchFaceLayers, // layers per face
nPatchPointLayers, // layers per point
@ -3050,8 +3066,8 @@ void Foam::autoLayerDriver::addLayers
addLayer.updateMesh
(
map,
identity(pp.size()),
identity(pp.nPoints())
identity(pp().size()),
identity(pp().nPoints())
);
// Collect layer faces and cells for outside loop.
@ -3098,7 +3114,7 @@ void Foam::autoLayerDriver::addLayers
(
addLayer,
meshQualityDict,
pp,
pp(),
newMesh,
patchDisp,
@ -3107,7 +3123,7 @@ void Foam::autoLayerDriver::addLayers
);
Info<< "Extruding " << countExtrusion(pp, extrudeStatus)
<< " out of " << returnReduce(pp.size(), sumOp<label>())
<< " out of " << returnReduce(pp().size(), sumOp<label>())
<< " faces. Removed extrusion at " << nTotChanged << " faces."
<< endl;
@ -3117,8 +3133,8 @@ void Foam::autoLayerDriver::addLayers
}
// Reset mesh points and start again
meshMover.movePoints(oldPoints);
meshMover.correct();
meshMover().movePoints(oldPoints);
meshMover().correct();
Info<< endl;
}
@ -3173,6 +3189,7 @@ void Foam::autoLayerDriver::addLayers
(
false,
false,
scalarField(mesh.nCells(), 1.0),
decomposer,
distributor
);
@ -3211,7 +3228,7 @@ void Foam::autoLayerDriver::doLayers
fvMeshDistribute& distributor
)
{
fvMesh& mesh = meshRefiner_.mesh();
const fvMesh& mesh = meshRefiner_.mesh();
Info<< nl
<< "Shrinking and layer addition phase" << nl
@ -3245,59 +3262,80 @@ void Foam::autoLayerDriver::doLayers
}
else
{
autoPtr<indirectPrimitivePatch> ppPtr
// Check that outside of mesh is not multiply connected.
checkMeshManifold();
// Check initial mesh
Info<< "Checking initial mesh ..." << endl;
labelHashSet wrongFaces(mesh.nFaces()/100);
motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
const label nInitErrors = returnReduce
(
meshRefinement::makePatch
(
mesh,
patchIDs
)
wrongFaces.size(),
sumOp<label>()
);
indirectPrimitivePatch& pp = ppPtr();
// Construct iterative mesh mover.
Info<< "Constructing mesh displacer ..." << endl;
Info<< "Detected " << nInitErrors << " illegal faces"
<< " (concave, zero area or negative cell pyramid volume)"
<< endl;
// Balance
if (Pstream::parRun())
{
const pointMesh& pMesh = pointMesh::New(mesh);
motionSmoother meshMover
(
mesh,
pp,
patchIDs,
meshRefinement::makeDisplacementField(pMesh, patchIDs),
motionDict
);
// Check that outside of mesh is not multiply connected.
checkMeshManifold();
// Check initial mesh
Info<< "Checking initial mesh ..." << endl;
labelHashSet wrongFaces(mesh.nFaces()/100);
motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
const label nInitErrors = returnReduce
(
wrongFaces.size(),
sumOp<label>()
);
Info<< "Detected " << nInitErrors << " illegal faces"
<< " (concave, zero area or negative cell pyramid volume)"
Info<< nl
<< "Doing initial balancing" << nl
<< "-----------------------" << nl
<< endl;
// Do all topo changes
addLayers
scalarField cellWeights(mesh.nCells(), 1);
forAll(numLayers, patchI)
{
if (numLayers[patchI] > 0)
{
const polyPatch& pp = mesh.boundaryMesh()[patchI];
forAll(pp.faceCells(), i)
{
cellWeights[pp.faceCells()[i]] += numLayers[patchI];
}
}
}
// Balance mesh (and meshRefinement). No restriction on face zones
// and baffles.
autoPtr<mapDistributePolyMesh> map = meshRefiner_.balance
(
layerParams,
motionDict,
nInitErrors,
meshMover,
false,
false,
cellWeights,
decomposer,
distributor
);
//{
// globalIndex globalCells(mesh.nCells());
//
// Info<< "** Distribution after balancing:" << endl;
// for (label procI = 0; procI < Pstream::nProcs(); procI++)
// {
// Info<< " " << procI << '\t'
// << globalCells.localSize(procI) << endl;
// }
// Info<< endl;
//}
}
// Do all topo changes
addLayers
(
layerParams,
motionDict,
patchIDs,
nInitErrors,
decomposer,
distributor
);
}
}

View File

@ -523,8 +523,8 @@ public:
(
const layerParameters& layerParams,
const dictionary& motionDict,
const labelList& patchIDs,
const label nAllowableErrors,
motionSmoother& meshMover,
decompositionMethod& decomposer,
fvMeshDistribute& distributor
);

View File

@ -35,7 +35,6 @@ License
#include "refinementSurfaces.H"
#include "shellSurfaces.H"
#include "mapDistributePolyMesh.H"
#include "mathConstants.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -210,7 +209,8 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
"feature refinement iteration " + name(iter),
decomposer_,
distributor_,
cellsToRefine
cellsToRefine,
refineParams.maxLoadUnbalance()
);
}
}
@ -310,7 +310,8 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
"surface refinement iteration " + name(iter),
decomposer_,
distributor_,
cellsToRefine
cellsToRefine,
refineParams.maxLoadUnbalance()
);
}
return iter;
@ -492,7 +493,8 @@ Foam::label Foam::autoRefineDriver::shellRefine
"shell refinement iteration " + name(iter),
decomposer_,
distributor_,
cellsToRefine
cellsToRefine,
refineParams.maxLoadUnbalance()
);
}
meshRefiner_.userFaceData().clear();
@ -776,11 +778,14 @@ void Foam::autoRefineDriver::doRefine
const_cast<Time&>(mesh.time())++;
}
// Do final balancing. Keep zoned faces on one processor.
// Do final balancing. Keep zoned faces on one processor since the
// snap phase will convert them to baffles and this only works for
// internal faces.
meshRefiner_.balance
(
true,
false,
scalarField(mesh.nCells(), 1), // dummy weights
decomposer_,
distributor_
);

View File

@ -43,7 +43,8 @@ Foam::refinementParameters::refinementParameters
minRefineCells_(readLabel(dict.lookup("minimumRefine"))),
curvature_(readScalar(dict.lookup("curvature"))),
nBufferLayers_(readLabel(dict.lookup("nBufferLayers"))),
keepPoints_(dict.lookup("keepPoints"))
keepPoints_(dict.lookup("keepPoints")),
maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance",0))
{}
@ -53,7 +54,8 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
maxLocalCells_(readLabel(dict.lookup("maxLocalCells"))),
minRefineCells_(readLabel(dict.lookup("minRefinementCells"))),
nBufferLayers_(readLabel(dict.lookup("nCellsBetweenLevels"))),
keepPoints_(pointField(1, dict.lookup("locationInMesh")))
keepPoints_(pointField(1, dict.lookup("locationInMesh"))),
maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance",0))
{
scalar featAngle(readScalar(dict.lookup("resolveFeatureAngle")));

View File

@ -73,6 +73,9 @@ class refinementParameters
//- Areas to keep
const pointField keepPoints_;
//- Allowed load unbalance
scalar maxLoadUnbalance_;
// Private Member Functions
@ -134,6 +137,12 @@ public:
return keepPoints_;
}
//- Allowed load unbalance
scalar maxLoadUnbalance() const
{
return maxLoadUnbalance_;
}
// Other

View File

@ -552,13 +552,16 @@ void Foam::meshRefinement::calcLocalRegions
const globalIndex& globalCells,
const labelList& globalRegion,
const Map<label>& coupledRegionToMaster,
const scalarField& cellWeights,
Map<label>& globalToLocalRegion,
pointField& localPoints
pointField& localPoints,
scalarField& localWeights
) const
{
globalToLocalRegion.resize(globalRegion.size());
DynamicList<point> localCc(globalRegion.size()/2);
DynamicList<scalar> localWts(globalRegion.size()/2);
forAll(globalRegion, cellI)
{
@ -573,6 +576,7 @@ void Foam::meshRefinement::calcLocalRegions
// I am master. Allocate region for me.
globalToLocalRegion.insert(globalRegion[cellI], localCc.size());
localCc.append(mesh_.cellCentres()[cellI]);
localWts.append(cellWeights[cellI]);
}
}
else
@ -581,11 +585,13 @@ void Foam::meshRefinement::calcLocalRegions
if (globalToLocalRegion.insert(globalRegion[cellI], localCc.size()))
{
localCc.append(mesh_.cellCentres()[cellI]);
localWts.append(cellWeights[cellI]);
}
}
}
localPoints.transfer(localCc);
localWeights.transfer(localWts);
if (localPoints.size() != globalToLocalRegion.size())
{
@ -924,6 +930,7 @@ Foam::label Foam::meshRefinement::countHits() const
// Determine distribution to move connected regions onto one processor.
Foam::labelList Foam::meshRefinement::decomposeCombineRegions
(
const scalarField& cellWeights,
const boolList& blockedFace,
const List<labelPair>& explicitConnections,
decompositionMethod& decomposer
@ -965,14 +972,17 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
Map<label> globalToLocalRegion;
pointField localPoints;
scalarField localWeights;
calcLocalRegions
(
globalCells,
globalRegion,
coupledRegionToMaster,
cellWeights,
globalToLocalRegion,
localPoints
localPoints,
localWeights
);
@ -984,7 +994,7 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
if (isA<geomDecomp>(decomposer))
{
regionDistribution = decomposer.decompose(localPoints);
regionDistribution = decomposer.decompose(localPoints, localWeights);
}
else
{
@ -998,7 +1008,12 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
regionRegions
);
regionDistribution = decomposer.decompose(regionRegions, localPoints);
regionDistribution = decomposer.decompose
(
regionRegions,
localPoints,
localWeights
);
}
@ -1058,6 +1073,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
(
const bool keepZoneFaces,
const bool keepBaffles,
const scalarField& cellWeights,
decompositionMethod& decomposer,
fvMeshDistribute& distributor
)
@ -1151,6 +1167,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
distribution = decomposeCombineRegions
(
cellWeights,
blockedFace,
couples,
decomposer
@ -1170,13 +1187,21 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
else
{
// Normal decomposition
distribution = decomposer.decompose(mesh_.cellCentres());
distribution = decomposer.decompose
(
mesh_.cellCentres(),
cellWeights
);
}
}
else
{
// Normal decomposition
distribution = decomposer.decompose(mesh_.cellCentres());
distribution = decomposer.decompose
(
mesh_.cellCentres(),
cellWeights
);
}
if (debug)

View File

@ -214,9 +214,11 @@ private:
const globalIndex& globalCells,
const labelList& globalRegion,
const Map<label>& coupledRegionToMaster,
const scalarField& cellWeights,
Map<label>& globalToLocalRegion,
pointField& localPoints
pointField& localPoints,
scalarField& localWeights
) const;
//- Convert region into global index.
@ -574,6 +576,7 @@ public:
// (e.g. to keep baffles together)
labelList decomposeCombineRegions
(
const scalarField& cellWeights,
const boolList& blockedFace,
const List<labelPair>& explicitConnections,
decompositionMethod&
@ -587,6 +590,7 @@ public:
(
const bool keepZoneFaces,
const bool keepBaffles,
const scalarField& cellWeights,
decompositionMethod& decomposer,
fvMeshDistribute& distributor
);
@ -648,7 +652,8 @@ public:
const string& msg,
decompositionMethod& decomposer,
fvMeshDistribute& distributor,
const labelList& cellsToRefine
const labelList& cellsToRefine,
const scalar maxLoadUnbalance
);

View File

@ -37,6 +37,7 @@ License
#include "mapDistributePolyMesh.H"
#include "featureEdgeMesh.H"
#include "Cloud.H"
#include "globalIndex.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -1244,40 +1245,122 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::refine
}
// Do refinement of consistent set of cells followed by truncation and
// load balancing.
//// Do refinement of consistent set of cells followed by truncation and
//// load balancing.
//Foam::autoPtr<Foam::mapDistributePolyMesh>
//Foam::meshRefinement::refineAndBalance
//(
// const string& msg,
// decompositionMethod& decomposer,
// fvMeshDistribute& distributor,
// const labelList& cellsToRefine
//)
//{
// // Do all refinement
// refine(cellsToRefine);
//
// if (debug)
// {
// Pout<< "Writing refined but unbalanced " << msg
// << " mesh to time " << timeName() << endl;
// write
// (
// debug,
// mesh_.time().path()
// /timeName()
// );
// Pout<< "Dumped debug data in = "
// << mesh_.time().cpuTimeIncrement() << " s" << endl;
//
// // test all is still synced across proc patches
// checkData();
// }
//
// Info<< "Refined mesh in = "
// << mesh_.time().cpuTimeIncrement() << " s" << endl;
// printMeshInfo(debug, "After refinement " + msg);
//
//
// // Load balancing
// // ~~~~~~~~~~~~~~
//
// autoPtr<mapDistributePolyMesh> distMap;
//
// if (Pstream::nProcs() > 1)
// {
// scalarField cellWeights(mesh_.nCells(), 1);
//
// distMap = balance
// (
// false, //keepZoneFaces
// false, //keepBaffles
// cellWeights,
// decomposer,
// distributor
// );
//
// Info<< "Balanced mesh in = "
// << mesh_.time().cpuTimeIncrement() << " s" << endl;
//
// printMeshInfo(debug, "After balancing " + msg);
//
//
// if (debug)
// {
// Pout<< "Writing balanced " << msg
// << " mesh to time " << timeName() << endl;
// write
// (
// debug,
// mesh_.time().path()/timeName()
// );
// Pout<< "Dumped debug data in = "
// << mesh_.time().cpuTimeIncrement() << " s" << endl;
//
// // test all is still synced across proc patches
// checkData();
// }
// }
//
// return distMap;
//}
// Do load balancing followed by refinement of consistent set of cells.
Foam::autoPtr<Foam::mapDistributePolyMesh>
Foam::meshRefinement::refineAndBalance
(
const string& msg,
decompositionMethod& decomposer,
fvMeshDistribute& distributor,
const labelList& cellsToRefine
const labelList& initCellsToRefine,
const scalar maxLoadUnbalance
)
{
// Do all refinement
refine(cellsToRefine);
labelList cellsToRefine(initCellsToRefine);
if (debug)
{
Pout<< "Writing refined but unbalanced " << msg
<< " mesh to time " << timeName() << endl;
write
(
debug,
mesh_.time().path()
/timeName()
);
Pout<< "Dumped debug data in = "
<< mesh_.time().cpuTimeIncrement() << " s" << endl;
// test all is still synced across proc patches
checkData();
}
Info<< "Refined mesh in = "
<< mesh_.time().cpuTimeIncrement() << " s" << endl;
printMeshInfo(debug, "After refinement " + msg);
//{
// globalIndex globalCells(mesh_.nCells());
//
// Info<< "** Distribution before balancing/refining:" << endl;
// for (label procI = 0; procI < Pstream::nProcs(); procI++)
// {
// Info<< " " << procI << '\t'
// << globalCells.localSize(procI) << endl;
// }
// Info<< endl;
//}
//{
// globalIndex globalCells(cellsToRefine.size());
//
// Info<< "** Cells to be refined:" << endl;
// for (label procI = 0; procI < Pstream::nProcs(); procI++)
// {
// Info<< " " << procI << '\t'
// << globalCells.localSize(procI) << endl;
// }
// Info<< endl;
//}
// Load balancing
@ -1287,20 +1370,62 @@ Foam::meshRefinement::refineAndBalance
if (Pstream::nProcs() > 1)
{
distMap = balance
// First check if we need to balance at all. Precalculate number of
// cells after refinement and see what maximum difference is.
scalar nNewCells = scalar(mesh_.nCells() + 7*cellsToRefine.size());
scalar nIdealNewCells =
returnReduce(nNewCells, sumOp<scalar>())
/ Pstream::nProcs();
scalar unbalance = returnReduce
(
false, //keepZoneFaces
false, //keepBaffles
decomposer,
distributor
mag(1.0-nNewCells/nIdealNewCells),
maxOp<scalar>()
);
Info<< "Balanced mesh in = "
<< mesh_.time().cpuTimeIncrement() << " s" << endl;
if (unbalance <= maxLoadUnbalance)
{
Info<< "Skipping balancing since max unbalance " << unbalance
<< " in = "
<< mesh_.time().cpuTimeIncrement() << " s" << endl;
}
else
{
scalarField cellWeights(mesh_.nCells(), 1);
forAll(cellsToRefine, i)
{
cellWeights[cellsToRefine[i]] += 7;
}
distMap = balance
(
false, //keepZoneFaces
false, //keepBaffles
cellWeights,
decomposer,
distributor
);
// Update cells to refine
distMap().distributeCellIndices(cellsToRefine);
Info<< "Balanced mesh in = "
<< mesh_.time().cpuTimeIncrement() << " s" << endl;
}
//{
// globalIndex globalCells(mesh_.nCells());
//
// Info<< "** Distribution after balancing:" << endl;
// for (label procI = 0; procI < Pstream::nProcs(); procI++)
// {
// Info<< " " << procI << '\t'
// << globalCells.localSize(procI) << endl;
// }
// Info<< endl;
//}
printMeshInfo(debug, "After balancing " + msg);
if (debug)
{
Pout<< "Writing balanced " << msg
@ -1318,6 +1443,46 @@ Foam::meshRefinement::refineAndBalance
}
}
// Refinement
// ~~~~~~~~~~
refine(cellsToRefine);
if (debug)
{
Pout<< "Writing refined " << msg
<< " mesh to time " << timeName() << endl;
write
(
debug,
mesh_.time().path()
/timeName()
);
Pout<< "Dumped debug data in = "
<< mesh_.time().cpuTimeIncrement() << " s" << endl;
// test all is still synced across proc patches
checkData();
}
Info<< "Refined mesh in = "
<< mesh_.time().cpuTimeIncrement() << " s" << endl;
//{
// globalIndex globalCells(mesh_.nCells());
//
// Info<< "** After refinement distribution:" << endl;
// for (label procI = 0; procI < Pstream::nProcs(); procI++)
// {
// Info<< " " << procI << '\t'
// << globalCells.localSize(procI) << endl;
// }
// Info<< endl;
//}
printMeshInfo(debug, "After refinement " + msg);
return distMap;
}