Skip to content
GitLab
Commit 01b785f2 authored by mattijs's avatar mattijs
Browse files

ENH: 2:1 constraint

parent 47bd5134
Hide whitespace changes
Inline Side-by-side
src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
View file @ 01b785f2
......@@ -1564,6 +1564,7 @@ void Foam::hexRef8::walkFaceFromMid
Foam::label Foam::hexRef8::faceConsistentRefinement
(
const bool maxSet,
const labelUList& cellLevel,
bitSet& refineCell
) const
{
......@@ -1573,10 +1574,10 @@ Foam::label Foam::hexRef8::faceConsistentRefinement
for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
{
label own = mesh_.faceOwner()[facei];
label nei = mesh_.faceNeighbour()[facei];
label ownLevel = cellLevel[own] + refineCell.get(own);
label ownLevel = cellLevel_[own] + refineCell.get(own);
label neiLevel = cellLevel_[nei] + refineCell.get(nei);
label nei = mesh_.faceNeighbour()[facei];
label neiLevel = cellLevel[nei] + refineCell.get(nei);
if (ownLevel > (neiLevel+1))
{
......@@ -1613,7 +1614,7 @@ Foam::label Foam::hexRef8::faceConsistentRefinement
{
label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
neiLevel[i] = cellLevel_[own] + refineCell.get(own);
neiLevel[i] = cellLevel[own] + refineCell.get(own);
}
// Swap to neighbour
......@@ -1623,7 +1624,7 @@ Foam::label Foam::hexRef8::faceConsistentRefinement
forAll(neiLevel, i)
{
label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
label ownLevel = cellLevel_[own] + refineCell.get(own);
label ownLevel = cellLevel[own] + refineCell.get(own);
if (ownLevel > (neiLevel[i]+1))
{
......@@ -1650,6 +1651,7 @@ Foam::label Foam::hexRef8::faceConsistentRefinement
// Debug: check if wanted refinement is compatible with 2:1
void Foam::hexRef8::checkWantedRefinementLevels
(
const labelUList& cellLevel,
const labelList& cellsToRefine
) const
{
......@@ -1658,10 +1660,10 @@ void Foam::hexRef8::checkWantedRefinementLevels
for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
{
label own = mesh_.faceOwner()[facei];
label nei = mesh_.faceNeighbour()[facei];
label ownLevel = cellLevel[own] + refineCell.get(own);
label ownLevel = cellLevel_[own] + refineCell.get(own);
label neiLevel = cellLevel_[nei] + refineCell.get(nei);
label nei = mesh_.faceNeighbour()[facei];
label neiLevel = cellLevel[nei] + refineCell.get(nei);
if (mag(ownLevel-neiLevel) > 1)
{
......@@ -1669,11 +1671,11 @@ void Foam::hexRef8::checkWantedRefinementLevels
dumpCell(nei);
FatalErrorInFunction
<< "cell:" << own
<< " current level:" << cellLevel_[own]
<< " current level:" << cellLevel[own]
<< " level after refinement:" << ownLevel
<< nl
<< "neighbour cell:" << nei
<< " current level:" << cellLevel_[nei]
<< " current level:" << cellLevel[nei]
<< " level after refinement:" << neiLevel
<< nl
<< "which does not satisfy 2:1 constraints anymore."
......@@ -1689,7 +1691,7 @@ void Foam::hexRef8::checkWantedRefinementLevels
{
label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
neiLevel[i] = cellLevel_[own] + refineCell.get(own);
neiLevel[i] = cellLevel[own] + refineCell.get(own);
}
// Swap to neighbour
......@@ -1701,7 +1703,7 @@ void Foam::hexRef8::checkWantedRefinementLevels
label facei = i + mesh_.nInternalFaces();
label own = mesh_.faceOwner()[facei];
label ownLevel = cellLevel_[own] + refineCell.get(own);
label ownLevel = cellLevel[own] + refineCell.get(own);
if (mag(ownLevel - neiLevel[i]) > 1)
{
......@@ -1715,7 +1717,7 @@ void Foam::hexRef8::checkWantedRefinementLevels
<< " on patch " << patchi << " "
<< mesh_.boundaryMesh()[patchi].name()
<< " owner cell " << own
<< " current level:" << cellLevel_[own]
<< " current level:" << cellLevel[own]
<< " level after refinement:" << ownLevel
<< nl
<< " (coupled) neighbour cell will get refinement "
......@@ -2251,6 +2253,7 @@ Foam::hexRef8::hexRef8
Foam::labelList Foam::hexRef8::consistentRefinement
(
const labelUList& cellLevel,
const labelList& cellsToRefine,
const bool maxSet
) const
......@@ -2264,7 +2267,12 @@ Foam::labelList Foam::hexRef8::consistentRefinement
while (true)
{
label nChanged = faceConsistentRefinement(maxSet, refineCell);
label nChanged = faceConsistentRefinement
(
maxSet,
cellLevel,
refineCell
);
reduce(nChanged, sumOp<label>());
......@@ -2286,7 +2294,7 @@ Foam::labelList Foam::hexRef8::consistentRefinement
if (debug)
{
checkWantedRefinementLevels(newCellsToRefine);
checkWantedRefinementLevels(cellLevel, newCellsToRefine);
}
return newCellsToRefine;
......@@ -3089,11 +3097,11 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
refineCell.set(celli);
}
}
faceConsistentRefinement(true, refineCell);
faceConsistentRefinement(true, cellLevel_, refineCell);
while (true)
{
label nChanged = faceConsistentRefinement(true, refineCell);
label nChanged = faceConsistentRefinement(true, cellLevel_, refineCell);
reduce(nChanged, sumOp<label>());
......@@ -3141,7 +3149,7 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
const bitSet savedRefineCell(refineCell);
label nChanged = faceConsistentRefinement(true, refineCell);
label nChanged = faceConsistentRefinement(true, cellLevel_, refineCell);
{
cellSet cellsOut2
......
src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.H
View file @ 01b785f2
......@@ -300,11 +300,16 @@ class hexRef8
label faceConsistentRefinement
(
const bool maxSet,
const labelUList& cellLevel,
bitSet& refineCell
) const;
//- Check wanted refinement for 2:1 consistency
void checkWantedRefinementLevels(const labelList&) const;
void checkWantedRefinementLevels
(
const labelUList& cellLevel,
const labelList&
) const;
// Cellshape recognition
......@@ -420,10 +425,25 @@ public:
// removes cells to refine (maxSet = false)
labelList consistentRefinement
(
const labelUList& cellLevel,
const labelList& cellsToRefine,
const bool maxSet
) const;
//- Given valid mesh and current cell level and proposed
// cells to refine calculate any clashes (due to 2:1) and return
// ok list of cells to refine.
// Either adds cells to refine to set (maxSet = true) or
// removes cells to refine (maxSet = false)
labelList consistentRefinement
(
const labelList& cellsToRefine,
const bool maxSet
) const
{
return consistentRefinement(cellLevel_, cellsToRefine, maxSet);
}
//- Like consistentRefinement but slower:
//
// - specify number of cells between consecutive refinement levels
......
src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
View file @ 01b785f2
......@@ -1037,6 +1037,16 @@ public:
const scalar maxLoadUnbalance
);
//- Calculate list of cells to directionally refine
labelList directionalRefineCandidates
(
const label maxGlobalCells,
const label maxLocalCells,
const labelList& currentLevel,
const direction dir
) const;
//- Directionally refine in direction cmpt
autoPtr<mapPolyMesh> directionalRefine
(
......
src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
View file @ 01b785f2
......@@ -2078,7 +2078,7 @@ Foam::labelList Foam::meshRefinement::refineCandidates
label nAllowRefine = labelMax / Pstream::nProcs();
// Marked for refinement (>= 0) or not (-1). Actual value is the
// index of the surface it intersects.
// index of the surface it intersects / shell it is inside.
labelList refineCell(mesh_.nCells(), -1);
label nRefine = 0;
......@@ -2629,6 +2629,106 @@ Foam::meshRefinement::balanceAndRefine
}
Foam::labelList Foam::meshRefinement::directionalRefineCandidates
(
const label maxGlobalCells,
const label maxLocalCells,
const labelList& currentLevel,
const direction dir
) const
{
const labelList& cellLevel = meshCutter_.cellLevel();
const pointField& cellCentres = mesh_.cellCentres();
label totNCells = mesh_.globalData().nTotalCells();
labelList cellsToRefine;
if (totNCells >= maxGlobalCells)
{
Info<< "No cells marked for refinement since reached limit "
<< maxGlobalCells << '.' << endl;
}
else
{
// Disable refinement shortcut. nAllowRefine is per processor limit.
label nAllowRefine = labelMax / Pstream::nProcs();
// Marked for refinement (>= 0) or not (-1). Actual value is the
// index of the surface it intersects / shell it is inside
labelList refineCell(mesh_.nCells(), -1);
label nRefine = 0;
// Find cells inside the shells with directional levels
labelList insideShell;
shells_.findDirectionalLevel
(
cellCentres,
cellLevel,
currentLevel, // current directional level
dir,
insideShell
);
// Mark for refinement
forAll(insideShell, celli)
{
if (insideShell[celli] >= 0)
{
bool reachedLimit = !markForRefine
(
insideShell[celli], // mark with any positive value
nAllowRefine,
refineCell[celli],
nRefine
);
if (reachedLimit)
{
if (debug)
{
Pout<< "Stopped refining cells"
<< " since reaching my cell limit of "
<< mesh_.nCells()+nAllowRefine << endl;
}
break;
}
}
}
// Limit refinement
// ~~~~~~~~~~~~~~~~
{
label nUnmarked = unmarkInternalRefinement(refineCell, nRefine);
if (nUnmarked > 0)
{
Info<< "Unmarked for refinement due to limit shells"
<< " : " << nUnmarked << " cells." << endl;
}
}
// Pack cells-to-refine
// ~~~~~~~~~~~~~~~~~~~~
cellsToRefine.setSize(nRefine);
nRefine = 0;
forAll(refineCell, cellI)
{
if (refineCell[cellI] != -1)
{
cellsToRefine[nRefine++] = cellI;
}
}
}
return cellsToRefine;
}
Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::directionalRefine
(
const string& msg,
......
src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
View file @ 01b785f2
......@@ -1572,190 +1572,6 @@ Foam::label Foam::snappyRefineDriver::shellRefine
return iter;
}
//XXXXXXXXXX
Foam::List<Foam::direction> Foam::snappyRefineDriver::faceDirection() const
{
const fvMesh& mesh = meshRefiner_.mesh();
List<direction> faceDir(mesh.nFaces());
vectorField nf(mesh.faceAreas());
nf /= mag(nf);
forAll(nf, facei)
{
const vector& n = nf[facei];
if (mag(n[0]) > 0.5)
{
faceDir[facei] = 0;
}
else if (mag(n[1]) > 0.5)
{
faceDir[facei] = 1;
}
else if (mag(n[2]) > 0.5)
{
faceDir[facei] = 2;
}
}
return faceDir;
}
Foam::label Foam::snappyRefineDriver::faceConsistentRefinement
(
const PackedBoolList& isDirFace,
const List<labelVector>& dirCellLevel,
const direction dir,
const bool maxSet,
PackedBoolList& refineCell
) const
{
const fvMesh& mesh = meshRefiner_.mesh();
label nChanged = 0;
// Internal faces.
for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
{
if (isDirFace[facei])
{
label own = mesh.faceOwner()[facei];
label ownLevel = dirCellLevel[own][dir] + refineCell.get(own);
label nei = mesh.faceNeighbour()[facei];
label neiLevel = dirCellLevel[nei][dir] + refineCell.get(nei);
if (ownLevel > (neiLevel+1))
{
if (maxSet)
{
refineCell.set(nei);
}
else
{
refineCell.unset(own);
}
nChanged++;
}
else if (neiLevel > (ownLevel+1))
{
if (maxSet)
{
refineCell.set(own);
}
else
{
refineCell.unset(nei);
}
nChanged++;
}
}
}
// Coupled faces. Swap owner level to get neighbouring cell level.
// (only boundary faces of neiLevel used)
labelList neiLevel(mesh.nFaces()-mesh.nInternalFaces());
forAll(neiLevel, i)
{
label own = mesh.faceOwner()[i+mesh.nInternalFaces()];
neiLevel[i] = dirCellLevel[own][dir] + refineCell.get(own);
}
// Swap to neighbour
syncTools::swapBoundaryFaceList(mesh, neiLevel);
// Now we have neighbour value see which cells need refinement
forAll(neiLevel, i)
{
label facei = i+mesh.nInternalFaces();
if (isDirFace[facei])
{
label own = mesh.faceOwner()[facei];
label ownLevel = dirCellLevel[own][dir] + refineCell.get(own);
if (ownLevel > (neiLevel[i]+1))
{
if (!maxSet)
{
refineCell.unset(own);
nChanged++;
}
}
else if (neiLevel[i] > (ownLevel+1))
{
if (maxSet)
{
refineCell.set(own);
nChanged++;
}
}
}
}
return nChanged;
}
Foam::labelList Foam::snappyRefineDriver::consistentDirectionalRefinement
(
const direction dir,
const List<labelVector>& dirCellLevel,
const labelUList& candidateCells,
const bool maxSet
) const
{
const fvMesh& mesh = meshRefiner_.mesh();
const List<direction> faceDir(faceDirection());
PackedBoolList isDirFace(mesh.nFaces());
forAll(faceDir, facei)
{
if (faceDir[facei] == dir)
{
isDirFace[facei] = true;
}
}
// Loop, modifying cellsToRefine, until no more changes to due to 2:1
// conflicts.
// maxSet = false : unselect cells to refine
// maxSet = true : select cells to refine
// Go to straight boolList.
PackedBoolList isRefineCell(mesh.nCells());
isRefineCell.set(candidateCells);
while (true)
{
label nChanged = faceConsistentRefinement
(
isDirFace,
dirCellLevel,
dir,
maxSet,
isRefineCell
);
reduce(nChanged, sumOp<label>());
if (debug)
{
Pout<< "snappyRefineDriver::consistentDirectionalRefinement"
<< " : Changed " << nChanged
<< " refinement levels due to 2:1 conflicts."
<< endl;
}
if (nChanged == 0)
{
break;
}
}
return isRefineCell.used();
}
Foam::label Foam::snappyRefineDriver::directionalShellRefine
......@@ -1814,53 +1630,31 @@ Foam::label Foam::snappyRefineDriver::directionalShellRefine
// - original cellLevel (using mapping) mentioned in levelIncrement
// - dirCellLevel not yet up to cellLevel+levelIncrement
labelList candidateCells;
// Extract component of directional level
labelList currentLevel(dirCellLevel.size());
forAll(dirCellLevel, celli)
{
// Extract component of directional level
labelList currentLevel(dirCellLevel.size());
forAll(dirCellLevel, celli)
{
currentLevel[celli] = dirCellLevel[celli][dir];
}
currentLevel[celli] = dirCellLevel[celli][dir];
}
// Find cells inside the shells
labelList insideShell;
shells.findDirectionalLevel
labelList candidateCells
(
meshRefiner_.directionalRefineCandidates
(
mesh.cellCentres(),
cellLevel,
currentLevel, // directional level
dir,
insideShell
);
// Extract
label n = 0;
forAll(insideShell, celli)
{
if (insideShell[celli] >= 0)
{
n++;
}
}
candidateCells.setSize(n);
n = 0;
forAll(insideShell, celli)
{
if (insideShell[celli] >= 0)
{
candidateCells[n++] = celli;
}
}
}
refineParams.maxGlobalCells(),
refineParams.maxLocalCells(),
currentLevel,
dir
)
);
// Extend to keep 2:1 ratio
labelList cellsToRefine
(
consistentDirectionalRefinement
meshRefiner_.meshCutter().consistentRefinement
(
dir,
dirCellLevel,
currentLevel,
candidateCells,
true
)
......