zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH update use of bitSet in dynamicRefineFvMesh #1075

Browse Source
This commit is contained in:
Mark Olesen
2018-11-16 10:30:02 +01:00
parent b4e26828d5
commit ad2baed5af
2 changed files with 206 additions and 295 deletions
Download Patch File Download Diff File Expand all files Collapse all files

471
src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
View File

@ -33,6 +33,7 @@ License
#include "pointFields.H" #include "pointFields.H"
#include "sigFpe.H" #include "sigFpe.H"
#include "cellSet.H" #include "cellSet.H"
#include "HashOps.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -44,32 +45,6 @@ namespace Foam
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
Foam::label Foam::dynamicRefineFvMesh::count
(
const bitSet& l,
const unsigned int val
)
{
label n = 0;
forAll(l, i)
{
if (l.get(i) == val)
{
n++;
}
// debug also serves to get-around Clang compiler trying to optimsie
// out this forAll loop under O3 optimisation
if (debug & 128)
{
Info<< "n=" << n << endl;
}
}
return n;
}
void Foam::dynamicRefineFvMesh::calculateProtectedCells void Foam::dynamicRefineFvMesh::calculateProtectedCells
( (
bitSet& unrefineableCell bitSet& unrefineableCell
@ -88,81 +63,89 @@ void Foam::dynamicRefineFvMesh::calculateProtectedCells
// Get neighbouring cell level // Get neighbouring cell level
labelList neiLevel(nFaces()-nInternalFaces()); labelList neiLevel(nFaces()-nInternalFaces());
for (label facei = nInternalFaces(); facei < nFaces(); facei++) for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
{ {
neiLevel[facei-nInternalFaces()] = cellLevel[faceOwner()[facei]]; neiLevel[facei-nInternalFaces()] = cellLevel[faceOwner()[facei]];
} }
syncTools::swapBoundaryFaceList(*this, neiLevel); syncTools::swapBoundaryFaceList(*this, neiLevel);
bitSet seedFace;
while (true) while (true)
{ {
// Pick up faces on border of protected cells // Pick up faces on border of protected cells
boolList seedFace(nFaces(), false); seedFace.reset();
seedFace.resize(nFaces());
forAll(faceNeighbour(), facei) for (label facei = 0; facei < nInternalFaces(); ++facei)
{ {
label own = faceOwner()[facei]; const label own = faceOwner()[facei];
bool ownProtected = unrefineableCell.get(own); const label nei = faceNeighbour()[facei];
label nei = faceNeighbour()[facei];
bool neiProtected = unrefineableCell.get(nei);
if (ownProtected && (cellLevel[nei] > cellLevel[own])) if
(
// Protected owner
(
unrefineableCell.test(own)
&& (cellLevel[nei] > cellLevel[own])
)
||
// Protected neighbour
(
unrefineableCell.test(nei)
&& (cellLevel[own] > cellLevel[nei])
)
)
{ {
seedFace[facei] = true; seedFace.set(facei);
}
else if (neiProtected && (cellLevel[own] > cellLevel[nei]))
{
seedFace[facei] = true;
} }
} }
for (label facei = nInternalFaces(); facei < nFaces(); facei++) for (label facei = nInternalFaces(); facei < nFaces(); facei++)
{ {
label own = faceOwner()[facei]; const label own = faceOwner()[facei];
bool ownProtected = unrefineableCell.get(own);
if if
( (
ownProtected // Protected owner
(
unrefineableCell.test(own)
&& (neiLevel[facei-nInternalFaces()] > cellLevel[own]) && (neiLevel[facei-nInternalFaces()] > cellLevel[own])
) )
)
{ {
seedFace[facei] = true; seedFace.set(facei);
} }
} }
syncTools::syncFaceList(*this, seedFace, orEqOp<bool>()); syncTools::syncFaceList(*this, seedFace, orEqOp<unsigned int>());
// Extend unrefineableCell // Extend unrefineableCell
bool hasExtended = false; bool hasExtended = false;
for (label facei = 0; facei < nInternalFaces(); facei++) for (label facei = 0; facei < nInternalFaces(); ++facei)
{ {
if (seedFace[facei]) if (seedFace.test(facei))
{ {
label own = faceOwner()[facei]; if (unrefineableCell.set(faceOwner()[facei]))
if (unrefineableCell.get(own) == 0)
{ {
unrefineableCell.set(own, 1);
hasExtended = true; hasExtended = true;
} }
if (unrefineableCell.set(faceNeighbour()[facei]))
{
hasExtended = true;
}
}
}
for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
{
if (seedFace.test(facei))
{
const label own = faceOwner()[facei];
label nei = faceNeighbour()[facei]; if (unrefineableCell.set(own))
if (unrefineableCell.get(nei) == 0)
{ {
unrefineableCell.set(nei, 1);
hasExtended = true;
}
}
}
for (label facei = nInternalFaces(); facei < nFaces(); facei++)
{
if (seedFace[facei])
{
label own = faceOwner()[facei];
if (unrefineableCell.get(own) == 0)
{
unrefineableCell.set(own, 1);
hasExtended = true; hasExtended = true;
} }
} }
@ -198,9 +181,9 @@ void Foam::dynamicRefineFvMesh::readDict()
// Rework into hashtable. // Rework into hashtable.
correctFluxes_.resize(fluxVelocities.size()); correctFluxes_.resize(fluxVelocities.size());
forAll(fluxVelocities, i) for (const auto& pr : fluxVelocities)
{ {
correctFluxes_.insert(fluxVelocities[i][0], fluxVelocities[i][1]); correctFluxes_.insert(pr.first(), pr.second());
} }
refineDict.readEntry("dumpLevel", dumpLevel_); refineDict.readEntry("dumpLevel", dumpLevel_);
@ -210,6 +193,7 @@ void Foam::dynamicRefineFvMesh::readDict()
void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm) void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm)
{ {
dynamicFvMesh::mapFields(mpm); dynamicFvMesh::mapFields(mpm);
// Correct the flux for modified/added faces. All the faces which only // Correct the flux for modified/added faces. All the faces which only
// have been renumbered will already have been handled by the mapping. // have been renumbered will already have been handled by the mapping.
{ {
@ -220,46 +204,42 @@ void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm)
// on the coarse cell that get split into four (or rather the // on the coarse cell that get split into four (or rather the
// master face gets modified and three faces get added from the master) // master face gets modified and three faces get added from the master)
// Estimate number of faces created // Estimate number of faces created
labelHashSet masterFaces
( bitSet masterFaces(nFaces());
max
(
mag(nFaces()-mpm.nOldFaces())/4,
nFaces()/100
)
);
forAll(faceMap, facei) forAll(faceMap, facei)
{ {
label oldFacei = faceMap[facei]; const label oldFacei = faceMap[facei];
if (oldFacei >= 0) if (oldFacei >= 0)
{ {
label masterFacei = reverseFaceMap[oldFacei]; const label masterFacei = reverseFaceMap[oldFacei];
if (masterFacei < 0) if (masterFacei < 0)
{ {
FatalErrorInFunction FatalErrorInFunction
<< "Problem: should not have removed faces" << "Problem: should not have removed faces"
<< " when refining." << " when refining."
<< nl << "face:" << facei << abort(FatalError); << nl << "face:" << facei << endl
<< abort(FatalError);
} }
else if (masterFacei != facei) else if (masterFacei != facei)
{ {
masterFaces.insert(masterFacei); masterFaces.set(masterFacei);
} }
} }
} }
if (debug) if (debug)
{ {
Pout<< "Found " << masterFaces.size() << " split faces " << endl; Pout<< "Found " << masterFaces.count() << " split faces " << endl;
} }
HashTable<surfaceScalarField*> fluxes HashTable<surfaceScalarField*> fluxes
( (
lookupClass<surfaceScalarField>() lookupClass<surfaceScalarField>()
); );
forAllIter(HashTable<surfaceScalarField*>, fluxes, iter) forAllIters(fluxes, iter)
{ {
if (!correctFluxes_.found(iter.key())) if (!correctFluxes_.found(iter.key()))
{ {
@ -282,13 +262,13 @@ void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm)
continue; continue;
} }
surfaceScalarField& phi = *iter();
if (UName == "NaN") if (UName == "NaN")
{ {
Pout<< "Setting surfaceScalarField " << iter.key() Pout<< "Setting surfaceScalarField " << iter.key()
<< " to NaN" << endl; << " to NaN" << endl;
surfaceScalarField& phi = *iter();
sigFpe::fillNan(phi.primitiveFieldRef()); sigFpe::fillNan(phi.primitiveFieldRef());
continue; continue;
@ -301,7 +281,6 @@ void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm)
<< endl; << endl;
} }
surfaceScalarField& phi = *iter();
const surfaceScalarField phiU const surfaceScalarField phiU
( (
fvc::interpolate fvc::interpolate
@ -312,9 +291,9 @@ void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm)
); );
// Recalculate new internal faces. // Recalculate new internal faces.
for (label facei = 0; facei < nInternalFaces(); facei++) for (label facei = 0; facei < nInternalFaces(); ++facei)
{ {
label oldFacei = faceMap[facei]; const label oldFacei = faceMap[facei];
if (oldFacei == -1) if (oldFacei == -1)
{ {
@ -329,8 +308,8 @@ void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm)
} }
// Recalculate new boundary faces. // Recalculate new boundary faces.
surfaceScalarField::Boundary& phiBf = surfaceScalarField::Boundary& phiBf = phi.boundaryFieldRef();
phi.boundaryFieldRef();
forAll(phiBf, patchi) forAll(phiBf, patchi)
{ {
fvsPatchScalarField& patchPhi = phiBf[patchi]; fvsPatchScalarField& patchPhi = phiBf[patchi];
@ -341,7 +320,7 @@ void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm)
forAll(patchPhi, i) forAll(patchPhi, i)
{ {
label oldFacei = faceMap[facei]; const label oldFacei = faceMap[facei];
if (oldFacei == -1) if (oldFacei == -1)
{ {
@ -354,23 +333,21 @@ void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm)
patchPhi[i] = patchPhiU[i]; patchPhi[i] = patchPhiU[i];
} }
facei++; ++facei;
} }
} }
// Update master faces // Update master faces
forAllConstIter(labelHashSet, masterFaces, iter) for (const label facei : masterFaces)
{ {
label facei = iter.key();
if (isInternalFace(facei)) if (isInternalFace(facei))
{ {
phi[facei] = phiU[facei]; phi[facei] = phiU[facei];
} }
else else
{ {
label patchi = boundaryMesh().whichPatch(facei); const label patchi = boundaryMesh().whichPatch(facei);
label i = facei - boundaryMesh()[patchi].start(); const label i = facei - boundaryMesh()[patchi].start();
const fvsPatchScalarField& patchPhiU = const fvsPatchScalarField& patchPhiU =
phiU.boundaryField()[patchi]; phiU.boundaryField()[patchi];
@ -424,9 +401,9 @@ Foam::dynamicRefineFvMesh::refine
if (debug) if (debug)
{ {
// Check map. // Check map.
for (label facei = 0; facei < nInternalFaces(); facei++) for (label facei = 0; facei < nInternalFaces(); ++facei)
{ {
label oldFacei = map().faceMap()[facei]; const label oldFacei = map().faceMap()[facei];
if (oldFacei >= nInternalFaces()) if (oldFacei >= nInternalFaces())
{ {
@ -474,14 +451,17 @@ Foam::dynamicRefineFvMesh::refine
forAll(newProtectedCell, celli) forAll(newProtectedCell, celli)
{ {
label oldCelli = map().cellMap()[celli]; const label oldCelli = map().cellMap()[celli];
newProtectedCell.set(celli, protectedCell_.get(oldCelli)); if (protectedCell_.test(oldCelli))
{
newProtectedCell.set(celli);
}
} }
protectedCell_.transfer(newProtectedCell); protectedCell_.transfer(newProtectedCell);
} }
// Debug: Check refinement levels (across faces only) // Debug: Check refinement levels (across faces only)
meshCutter_.checkRefinementLevels(-1, labelList(0)); meshCutter_.checkRefinementLevels(-1, labelList());
return map; return map;
} }
@ -508,21 +488,19 @@ Foam::dynamicRefineFvMesh::unrefine
Map<label> faceToSplitPoint(3*splitPoints.size()); Map<label> faceToSplitPoint(3*splitPoints.size());
{ {
forAll(splitPoints, i) for (const label pointi : splitPoints)
{ {
label pointi = splitPoints[i];
const labelList& pEdges = pointEdges()[pointi]; const labelList& pEdges = pointEdges()[pointi];
forAll(pEdges, j) for (const label edgei : pEdges)
{ {
label otherPointi = edges()[pEdges[j]].otherVertex(pointi); const label otherPointi = edges()[edgei].otherVertex(pointi);
const labelList& pFaces = pointFaces()[otherPointi]; const labelList& pFaces = pointFaces()[otherPointi];
forAll(pFaces, pFacei) for (const label facei : pFaces)
{ {
faceToSplitPoint.insert(pFaces[pFacei], otherPointi); faceToSplitPoint.insert(facei, otherPointi);
} }
} }
} }
@ -565,7 +543,7 @@ Foam::dynamicRefineFvMesh::unrefine
( (
lookupClass<surfaceScalarField>() lookupClass<surfaceScalarField>()
); );
forAllIter(HashTable<surfaceScalarField*>, fluxes, iter) forAllIters(fluxes, iter)
{ {
if (!correctFluxes_.found(iter.key())) if (!correctFluxes_.found(iter.key()))
{ {
@ -588,12 +566,11 @@ Foam::dynamicRefineFvMesh::unrefine
continue; continue;
} }
if (debug) DebugInfo
{ << "Mapping flux " << iter.key()
Info<< "Mapping flux " << iter.key()
<< " using interpolated flux " << UName << " using interpolated flux " << UName
<< endl; << endl;
}
surfaceScalarField& phi = *iter(); surfaceScalarField& phi = *iter();
surfaceScalarField::Boundary& phiBf = surfaceScalarField::Boundary& phiBf =
@ -609,15 +586,15 @@ Foam::dynamicRefineFvMesh::unrefine
); );
forAllConstIter(Map<label>, faceToSplitPoint, iter) forAllConstIters(faceToSplitPoint, iter)
{ {
label oldFacei = iter.key(); const label oldFacei = iter.key();
label oldPointi = iter(); const label oldPointi = iter.object();
if (reversePointMap[oldPointi] < 0) if (reversePointMap[oldPointi] < 0)
{ {
// midpoint was removed. See if face still exists. // midpoint was removed. See if face still exists.
label facei = reverseFaceMap[oldFacei]; const label facei = reverseFaceMap[oldFacei];
if (facei >= 0) if (facei >= 0)
{ {
@ -652,17 +629,17 @@ Foam::dynamicRefineFvMesh::unrefine
forAll(newProtectedCell, celli) forAll(newProtectedCell, celli)
{ {
label oldCelli = map().cellMap()[celli]; const label oldCelli = map().cellMap()[celli];
if (oldCelli >= 0) if (protectedCell_.test(oldCelli))
{ {
newProtectedCell.set(celli, protectedCell_.get(oldCelli)); newProtectedCell.set(celli);
} }
} }
protectedCell_.transfer(newProtectedCell); protectedCell_.transfer(newProtectedCell);
} }
// Debug: Check refinement levels (across faces only) // Debug: Check refinement levels (across faces only)
meshCutter_.checkRefinementLevels(-1, labelList(0)); meshCutter_.checkRefinementLevels(-1, labelList());
return map; return map;
} }
@ -677,9 +654,9 @@ Foam::dynamicRefineFvMesh::maxPointField(const scalarField& pFld) const
{ {
const labelList& pCells = pointCells()[pointi]; const labelList& pCells = pointCells()[pointi];
forAll(pCells, i) for (const label celli : pCells)
{ {
vFld[pCells[i]] = max(vFld[pCells[i]], pFld[pointi]); vFld[celli] = max(vFld[celli], pFld[pointi]);
} }
} }
return vFld; return vFld;
@ -695,9 +672,9 @@ Foam::dynamicRefineFvMesh::maxCellField(const volScalarField& vFld) const
{ {
const labelList& pCells = pointCells()[pointi]; const labelList& pCells = pointCells()[pointi];
forAll(pCells, i) for (const label celli : pCells)
{ {
pFld[pointi] = max(pFld[pointi], vFld[pCells[i]]); pFld[pointi] = max(pFld[pointi], vFld[celli]);
} }
} }
return pFld; return pFld;
@ -714,9 +691,9 @@ Foam::dynamicRefineFvMesh::cellToPoint(const scalarField& vFld) const
const labelList& pCells = pointCells()[pointi]; const labelList& pCells = pointCells()[pointi];
scalar sum = 0.0; scalar sum = 0.0;
forAll(pCells, i) for (const label celli : pCells)
{ {
sum += vFld[pCells[i]]; sum += vFld[celli];
} }
pFld[pointi] = sum/pCells.size(); pFld[pointi] = sum/pCells.size();
} }
@ -731,7 +708,7 @@ Foam::scalarField Foam::dynamicRefineFvMesh::error
const scalar maxLevel const scalar maxLevel
) const ) const
{ {
scalarField c(fld.size(), -1); scalarField c(fld.size(), scalar(-1));
forAll(fld, i) forAll(fld, i)
{ {
@ -774,7 +751,7 @@ void Foam::dynamicRefineFvMesh::selectRefineCandidates
{ {
if (cellError[celli] > 0) if (cellError[celli] > 0)
{ {
candidateCell.set(celli, 1); candidateCell.set(celli);
} }
} }
} }
@ -798,7 +775,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
calculateProtectedCells(unrefineableCell); calculateProtectedCells(unrefineableCell);
// Count current selection // Count current selection
label nLocalCandidates = count(candidateCell, 1); label nLocalCandidates = candidateCell.count();
label nCandidates = returnReduce(nLocalCandidates, sumOp<label>()); label nCandidates = returnReduce(nLocalCandidates, sumOp<label>());
// Collect all cells // Collect all cells
@ -806,16 +783,12 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
if (nCandidates < nTotToRefine) if (nCandidates < nTotToRefine)
{ {
forAll(candidateCell, celli) for (const label celli : candidateCell)
{ {
if if
( (
cellLevel[celli] < maxRefinement (!unrefineableCell.test(celli))
&& candidateCell.get(celli) && cellLevel[celli] < maxRefinement
&& (
unrefineableCell.empty()
|| !unrefineableCell.get(celli)
)
) )
{ {
candidates.append(celli); candidates.append(celli);
@ -825,18 +798,14 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
else else
{ {
// Sort by error? For now just truncate. // Sort by error? For now just truncate.
for (label level = 0; level < maxRefinement; level++) for (label level = 0; level < maxRefinement; ++level)
{ {
forAll(candidateCell, celli) for (const label celli : candidateCell)
{ {
if if
( (
cellLevel[celli] == level (!unrefineableCell.test(celli))
&& candidateCell.get(celli) && cellLevel[celli] == level
&& (
unrefineableCell.empty()
|| !unrefineableCell.get(celli)
)
) )
{ {
candidates.append(celli); candidates.append(celli);
@ -889,17 +858,13 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
if (protectedCell_.size()) if (protectedCell_.size())
{ {
// Get all points on a protected cell // Get all points on a protected cell
forAll(pointCells, pointI) forAll(pointCells, pointi)
{ {
const labelList& pCells = pointCells[pointI]; for (const label celli : pointCells[pointi])
forAll(pCells, pCellI)
{ {
label cellI = pCells[pCellI]; if (protectedCell_.test(celli))
if (protectedCell_[cellI])
{ {
protectedPoint[pointI] = true; protectedPoint.set(pointi);
break; break;
} }
} }
@ -910,36 +875,29 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
*this, *this,
protectedPoint, protectedPoint,
orEqOp<unsigned int>(), orEqOp<unsigned int>(),
0U 0u
); );
if (debug) DebugInfo<< "From "
{
Info<< "From "
<< returnReduce(protectedCell_.count(), sumOp<label>()) << returnReduce(protectedCell_.count(), sumOp<label>())
<< " protected cells found " << " protected cells found "
<< returnReduce(protectedPoint.count(), sumOp<label>()) << returnReduce(protectedPoint.count(), sumOp<label>())
<< " protected points." << endl; << " protected points." << endl;
} }
}
DynamicList<label> newSplitPoints(splitPoints.size()); DynamicList<label> newSplitPoints(splitPoints.size());
forAll(splitPoints, i) for (const label pointi : splitPoints)
{ {
label pointi = splitPoints[i];
if (!protectedPoint[pointi] && pFld[pointi] < unrefineLevel) if (!protectedPoint[pointi] && pFld[pointi] < unrefineLevel)
{ {
// Check that all cells are not marked // Check that all cells are not marked
const labelList& pCells = pointCells[pointi];
bool hasMarked = false; bool hasMarked = false;
forAll(pCells, pCelli) for (const label celli : pointCells[pointi])
{ {
if (markedCell.get(pCells[pCelli])) if (markedCell.test(celli))
{ {
hasMarked = true; hasMarked = true;
break; break;
@ -980,37 +938,29 @@ void Foam::dynamicRefineFvMesh::extendMarkedCells
) const ) const
{ {
// Mark faces using any marked cell // Mark faces using any marked cell
boolList markedFace(nFaces(), false); bitSet markedFace(nFaces());
forAll(markedCell, celli) for (const label celli : markedCell)
{ {
if (markedCell.get(celli)) markedFace.set(cells()[celli]); // set multiple faces
{
const cell& cFaces = cells()[celli];
forAll(cFaces, i)
{
markedFace[cFaces[i]] = true;
}
}
} }
syncTools::syncFaceList(*this, markedFace, orEqOp<bool>()); syncTools::syncFaceList(*this, markedFace, orEqOp<unsigned int>());
// Update cells using any markedFace // Update cells using any markedFace
for (label facei = 0; facei < nInternalFaces(); facei++) for (label facei = 0; facei < nInternalFaces(); ++facei)
{ {
if (markedFace[facei]) if (markedFace.test(facei))
{ {
markedCell.set(faceOwner()[facei], 1); markedCell.set(faceOwner()[facei]);
markedCell.set(faceNeighbour()[facei], 1); markedCell.set(faceNeighbour()[facei]);
} }
} }
for (label facei = nInternalFaces(); facei < nFaces(); facei++) for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
{ {
if (markedFace[facei]) if (markedFace.test(facei))
{ {
markedCell.set(faceOwner()[facei], 1); markedCell.set(faceOwner()[facei]);
} }
} }
} }
@ -1018,37 +968,31 @@ void Foam::dynamicRefineFvMesh::extendMarkedCells
void Foam::dynamicRefineFvMesh::checkEightAnchorPoints void Foam::dynamicRefineFvMesh::checkEightAnchorPoints
( (
bitSet& protectedCell, bitSet& protectedCell
label& nProtected
) const ) const
{ {
const labelList& cellLevel = meshCutter_.cellLevel(); const labelList& cellLevel = meshCutter_.cellLevel();
const labelList& pointLevel = meshCutter_.pointLevel(); const labelList& pointLevel = meshCutter_.pointLevel();
labelList nAnchorPoints(nCells(), 0); labelList nAnchorPoints(nCells(), Zero);
forAll(pointLevel, pointi) forAll(pointLevel, pointi)
{ {
const labelList& pCells = pointCells(pointi); const labelList& pCells = pointCells(pointi);
forAll(pCells, pCelli) for (const label celli : pCells)
{ {
label celli = pCells[pCelli];
if (pointLevel[pointi] <= cellLevel[celli]) if (pointLevel[pointi] <= cellLevel[celli])
{ {
// Check if cell has already 8 anchor points -> protect cell // Check if cell has already 8 anchor points -> protect cell
if (nAnchorPoints[celli] == 8) if (nAnchorPoints[celli] == 8)
{ {
if (protectedCell.set(celli, true)) protectedCell.set(celli);
{
nProtected++;
}
} }
if (!protectedCell[celli]) if (!protectedCell.test(celli))
{ {
nAnchorPoints[celli]++; ++nAnchorPoints[celli];
} }
} }
} }
@ -1057,10 +1001,9 @@ void Foam::dynamicRefineFvMesh::checkEightAnchorPoints
forAll(protectedCell, celli) forAll(protectedCell, celli)
{ {
if (!protectedCell[celli] && nAnchorPoints[celli] != 8) if (nAnchorPoints[celli] != 8)
{ {
protectedCell.set(celli, true); protectedCell.set(celli);
nProtected++;
} }
} }
} }
@ -1072,9 +1015,9 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
: :
dynamicFvMesh(io), dynamicFvMesh(io),
meshCutter_(*this), meshCutter_(*this),
dumpLevel_(false), protectedCell_(nCells()),
nRefinementIterations_(0), nRefinementIterations_(0),
protectedCell_(nCells(), 0) dumpLevel_(false)
{ {
// Read static part of dictionary // Read static part of dictionary
readDict(); readDict();
@ -1091,28 +1034,23 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
// Count number of points <= cellLevel // Count number of points <= cellLevel
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
labelList nAnchors(nCells(), 0); labelList nAnchors(nCells(), Zero);
label nProtected = 0;
forAll(pointCells(), pointi) forAll(pointCells(), pointi)
{ {
const labelList& pCells = pointCells()[pointi]; const labelList& pCells = pointCells()[pointi];
forAll(pCells, i) for (const label celli : pCells)
{ {
label celli = pCells[i]; if (!protectedCell_.test(celli))
if (!protectedCell_.get(celli))
{ {
if (pointLevel[pointi] <= cellLevel[celli]) if (pointLevel[pointi] <= cellLevel[celli])
{ {
nAnchors[celli]++; ++nAnchors[celli];
if (nAnchors[celli] > 8) if (nAnchors[celli] > 8)
{ {
protectedCell_.set(celli, 1); protectedCell_.set(celli);
nProtected++;
} }
} }
} }
@ -1128,22 +1066,22 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
{ {
labelList neiLevel(nFaces()); labelList neiLevel(nFaces());
for (label facei = 0; facei < nInternalFaces(); facei++) for (label facei = 0; facei < nInternalFaces(); ++facei)
{ {
neiLevel[facei] = cellLevel[faceNeighbour()[facei]]; neiLevel[facei] = cellLevel[faceNeighbour()[facei]];
} }
for (label facei = nInternalFaces(); facei < nFaces(); facei++) for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
{ {
neiLevel[facei] = cellLevel[faceOwner()[facei]]; neiLevel[facei] = cellLevel[faceOwner()[facei]];
} }
syncTools::swapFaceList(*this, neiLevel); syncTools::swapFaceList(*this, neiLevel);
boolList protectedFace(nFaces(), false); bitSet protectedFace(nFaces());
forAll(faceOwner(), facei) forAll(faceOwner(), facei)
{ {
label faceLevel = max const label faceLevel = max
( (
cellLevel[faceOwner()[facei]], cellLevel[faceOwner()[facei]],
neiLevel[facei] neiLevel[facei]
@ -1153,39 +1091,36 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
label nAnchors = 0; label nAnchors = 0;
forAll(f, fp) for (const label pointi : f)
{ {
if (pointLevel[f[fp]] <= faceLevel) if (pointLevel[pointi] <= faceLevel)
{ {
nAnchors++; ++nAnchors;
if (nAnchors > 4) if (nAnchors > 4)
{ {
protectedFace[facei] = true; protectedFace.set(facei);
break; break;
} }
} }
} }
} }
syncTools::syncFaceList(*this, protectedFace, orEqOp<bool>()); syncTools::syncFaceList(*this, protectedFace, orEqOp<unsigned int>());
for (label facei = 0; facei < nInternalFaces(); facei++) for (label facei = 0; facei < nInternalFaces(); ++facei)
{ {
if (protectedFace[facei]) if (protectedFace.test(facei))
{ {
protectedCell_.set(faceOwner()[facei], 1); protectedCell_.set(faceOwner()[facei]);
nProtected++; protectedCell_.set(faceNeighbour()[facei]);
protectedCell_.set(faceNeighbour()[facei], 1);
nProtected++;
} }
} }
for (label facei = nInternalFaces(); facei < nFaces(); facei++) for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
{ {
if (protectedFace[facei]) if (protectedFace.test(facei))
{ {
protectedCell_.set(faceOwner()[facei], 1); protectedCell_.set(faceOwner()[facei]);
nProtected++;
} }
} }
@ -1196,21 +1131,15 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
if (cFaces.size() < 6) if (cFaces.size() < 6)
{ {
if (protectedCell_.set(celli, 1)) protectedCell_.set(celli);
{
nProtected++;
}
} }
else else
{ {
forAll(cFaces, cFacei) for (const label cfacei : cFaces)
{ {
if (faces()[cFaces[cFacei]].size() < 4) if (faces()[cfacei].size() < 4)
{ {
if (protectedCell_.set(celli, 1)) protectedCell_.set(celli);
{
nProtected++;
}
break; break;
} }
} }
@ -1218,26 +1147,24 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
} }
// Check cells for 8 corner points // Check cells for 8 corner points
checkEightAnchorPoints(protectedCell_, nProtected); checkEightAnchorPoints(protectedCell_);
} }
if (returnReduce(nProtected, sumOp<label>()) == 0) if (!returnReduce(protectedCell_.any(), orOp<bool>()))
{ {
protectedCell_.clear(); protectedCell_.clear();
} }
else else
{ {
cellSet protectedCells
(
*this,
"protectedCells",
HashSetOps::used(protectedCell_)
);
cellSet protectedCells(*this, "protectedCells", nProtected); Info<< "Detected "
forAll(protectedCell_, celli) << returnReduce(protectedCells.size(), sumOp<label>())
{
if (protectedCell_[celli])
{
protectedCells.insert(celli);
}
}
Info<< "Detected " << returnReduce(nProtected, sumOp<label>())
<< " cells that are protected from refinement." << " cells that are protected from refinement."
<< " Writing these to cellSet " << " Writing these to cellSet "
<< protectedCells.name() << protectedCells.name()
@ -1248,12 +1175,6 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
} }
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::dynamicRefineFvMesh::~dynamicRefineFvMesh()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::dynamicRefineFvMesh::update() bool Foam::dynamicRefineFvMesh::update()
@ -1277,7 +1198,7 @@ bool Foam::dynamicRefineFvMesh::update()
).optionalSubDict(typeName + "Coeffs") ).optionalSubDict(typeName + "Coeffs")
); );
label refineInterval = refineDict.get<label>("refineInterval"); const label refineInterval = refineDict.get<label>("refineInterval");
bool hasChanged = false; bool hasChanged = false;
@ -1297,14 +1218,12 @@ bool Foam::dynamicRefineFvMesh::update()
} }
// Note: cannot refine at time 0 since no V0 present since mesh not // Note: cannot refine at time 0 since no V0 present since mesh not
// moved yet. // moved yet.
if (time().timeIndex() > 0 && time().timeIndex() % refineInterval == 0) if (time().timeIndex() > 0 && time().timeIndex() % refineInterval == 0)
{ {
label maxCells = refineDict.get<label>("maxCells"); const label maxCells = refineDict.get<label>("maxCells");
if (maxCells <= 0) if (maxCells <= 0)
{ {
@ -1315,7 +1234,7 @@ bool Foam::dynamicRefineFvMesh::update()
<< exit(FatalError); << exit(FatalError);
} }
label maxRefinement = refineDict.get<label>("maxRefinement"); const label maxRefinement = refineDict.get<label>("maxRefinement");
if (maxRefinement <= 0) if (maxRefinement <= 0)
{ {
@ -1339,8 +1258,7 @@ bool Foam::dynamicRefineFvMesh::update()
"unrefineLevel", "unrefineLevel",
GREAT GREAT
); );
const label nBufferLayers = const label nBufferLayers = refineDict.get<label>("nBufferLayers");
refineDict.get<label>("nBufferLayers");
// Cells marked for refinement or otherwise protected from unrefinement. // Cells marked for refinement or otherwise protected from unrefinement.
bitSet refineCell(nCells()); bitSet refineCell(nCells());
@ -1368,7 +1286,7 @@ bool Foam::dynamicRefineFvMesh::update()
) )
); );
label nCellsToRefine = returnReduce const label nCellsToRefine = returnReduce
( (
cellsToRefine.size(), sumOp<label>() cellsToRefine.size(), sumOp<label>()
); );
@ -1388,19 +1306,16 @@ bool Foam::dynamicRefineFvMesh::update()
forAll(cellMap, celli) forAll(cellMap, celli)
{ {
label oldCelli = cellMap[celli]; const label oldCelli = cellMap[celli];
if (oldCelli < 0) if
(
(oldCelli < 0)
|| (reverseCellMap[oldCelli] != celli)
|| (refineCell.test(oldCelli))
)
{ {
newRefineCell.set(celli, 1); newRefineCell.set(celli);
}
else if (reverseCellMap[oldCelli] != celli)
{
newRefineCell.set(celli, 1);
}
else
{
newRefineCell.set(celli, refineCell.get(oldCelli));
} }
} }
refineCell.transfer(newRefineCell); refineCell.transfer(newRefineCell);
@ -1408,7 +1323,7 @@ bool Foam::dynamicRefineFvMesh::update()
// Extend with a buffer layer to prevent neighbouring points // Extend with a buffer layer to prevent neighbouring points
// being unrefined. // being unrefined.
for (label i = 0; i < nBufferLayers; i++) for (label i = 0; i < nBufferLayers; ++i)
{ {
extendMarkedCells(refineCell); extendMarkedCells(refineCell);
} }
@ -1430,7 +1345,7 @@ bool Foam::dynamicRefineFvMesh::update()
) )
); );
label nSplitPoints = returnReduce const label nSplitPoints = returnReduce
( (
pointsToUnrefine.size(), pointsToUnrefine.size(),
sumOp<label>() sumOp<label>()
@ -1498,7 +1413,7 @@ bool Foam::dynamicRefineFvMesh::writeObject
false false
), ),
*this, *this,
dimensionedScalar("level", dimless, 0) dimensionedScalar(dimless, Zero)
); );
const labelList& cellLevel = meshCutter_.cellLevel(); const labelList& cellLevel = meshCutter_.cellLevel();

28
src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.H
View File

@ -29,7 +29,8 @@ Description
Determines which cells to refine/unrefine and does all in update(). Determines which cells to refine/unrefine and does all in update().
\verbatim
{
// How often to refine // How often to refine
refineInterval 1; refineInterval 1;
// Field to be refinement on // Field to be refinement on
@ -59,6 +60,8 @@ Description
// Write the refinement level as a volScalarField // Write the refinement level as a volScalarField
dumpLevel true; dumpLevel true;
}
\endverbatim
SourceFiles SourceFiles
@ -91,24 +94,21 @@ protected:
//- Mesh cutting engine //- Mesh cutting engine
hexRef8 meshCutter_; hexRef8 meshCutter_;
//- Dump cellLevel for post-processing
bool dumpLevel_;
//- Fluxes to map //- Fluxes to map
HashTable<word> correctFluxes_; HashTable<word> correctFluxes_;
//- Number of refinement/unrefinement steps done so far.
label nRefinementIterations_;
//- Protected cells (usually since not hexes) //- Protected cells (usually since not hexes)
bitSet protectedCell_; bitSet protectedCell_;
//- Number of refinement/unrefinement steps done so far.
label nRefinementIterations_;
//- Dump cellLevel for post-processing
bool dumpLevel_;
// Protected Member Functions // Protected Member Functions
//- Count set/unset elements in packedlist.
static label count(const bitSet&, const unsigned int);
//- Calculate cells that cannot be refined since would trigger //- Calculate cells that cannot be refined since would trigger
// refinement of protectedCell_ (since 2:1 refinement cascade) // refinement of protectedCell_ (since 2:1 refinement cascade)
void calculateProtectedCells(bitSet& unrefineableCell) const; void calculateProtectedCells(bitSet& unrefineableCell) const;
@ -180,11 +180,7 @@ protected:
void extendMarkedCells(bitSet& markedCell) const; void extendMarkedCells(bitSet& markedCell) const;
//- Check all cells have 8 anchor points //- Check all cells have 8 anchor points
void checkEightAnchorPoints void checkEightAnchorPoints(bitSet& protectedCell) const;
(
bitSet& protectedCell,
label& nProtected
) const;
//- Map single non-flux surface<Type>Field //- Map single non-flux surface<Type>Field
// for new internal faces (e.g. AMR refine). This currently // for new internal faces (e.g. AMR refine). This currently
@ -233,7 +229,7 @@ public:
//- Destructor //- Destructor
virtual ~dynamicRefineFvMesh(); virtual ~dynamicRefineFvMesh() = default;
// Member Functions // Member Functions