Commit 909b5906 authored by mattijs's avatar mattijs
Browse files

balance before refinement

parent 2f9e3fa6
......@@ -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
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
autoPtr<indirectPrimitivePatch> pp
(
meshRefinement::makePatch
(
mesh,
patchIDs
)
);
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
// Construct iterative mesh mover.
Info<< "Constructing mesh displacer ..." << endl;
autoPtr<motionSmoother> meshMover
(
new motionSmoother
(
mesh.edges(),
mesh.pointEdges()[v0],
v0,
v1
);
}
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;
{
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();
Info<< "Detected " << nInitErrors << " illegal faces"
<< " (concave, zero area or negative cell pyramid volume)"
<< endl;
// 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)"
// Balance
if (Pstream::parRun())
{
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
);
}
}
......
......@@ -523,8 +523,8 @@ public:
(
const layerParameters& layerParams,
const dictionary& motionDict,
const labelList& patchIDs,
const label nAllowableErrors,
motionSmoother& meshMover,
decompositionMethod& decomposer,
fvMeshDistribute& distributor
);
......
......@@ -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_
);
......
......@@ -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")));
......
......@@ -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
......
......@@ -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
);