Update mapping of cell data on unrefinement to include all neighbours

This commit is contained in:
Hrvoje Jasak 2019-09-09 15:06:20 +01:00
parent 814f3de0a5
commit 326c61d8ab
10 changed files with 304 additions and 219 deletions

View file

@ -147,11 +147,15 @@ int main(int argc, char *argv[])
// Topo changes container
directTopoChange meshMod(mesh);
// Dummy point region master
labelList pointRegionMaster(cellRegionMaster.size(), -1);
// Insert mesh refinement into directTopoChange.
faceRemover.setRefinement
(
facesToRemove,
cellRegion,
pointRegionMaster,
cellRegionMaster,
meshMod
);

View file

@ -5345,10 +5345,17 @@ void Foam::hexRef8::setUnrefinement
// Insert all commands to combine cells. Never fails so don't have to
// test for success.
// Added dummy pointRegionMaster, which will be ignored.
// Probably, splitPointLabels should be used
// HJ, 6/Sep/2019
labelList pointRegionMaster(cellRegionMaster.size(), label(-1));
faceRemover_.setRefinement
(
facesToRemove,
cellRegion,
pointRegionMaster,
cellRegionMaster,
meshMod
);

View file

@ -149,10 +149,14 @@ Foam::Xfer<Foam::boolList> Foam::removeFaces::affectedFaces
{
const label& region = cellRegion[cellI];
if (region != -1 && (cellI != cellRegionMaster[region]))
// Bug fix: map coarse cell from all fine neighbours
// All cells will be removed and replacew with a new one from
// the central point. HJ, 9/Sep/2019
if (region != -1)
{
// Get this cell (list of cell faces) and mark all of its faces
const labelList& cFaces = meshCells[cellI];
forAll(cFaces, cFaceI)
{
affectedFace[cFaces[cFaceI]] = true;
@ -284,7 +288,7 @@ Foam::label Foam::removeFaces::compatibleRemoves
(
const labelList& facesToRemove,
labelList& cellRegion,
labelList& regionMaster,
labelList& cellRegionMaster,
labelList& newFacesToRemove
) const
{
@ -298,8 +302,8 @@ Foam::label Foam::removeFaces::compatibleRemoves
cellRegion.setSize(nCells);
cellRegion = -1;
regionMaster.setSize(nCells);
regionMaster = -1;
cellRegionMaster.setSize(nCells);
cellRegionMaster = -1;
// Count regions
label nRegions = 0;
@ -342,7 +346,7 @@ Foam::label Foam::removeFaces::compatibleRemoves
// Make owner (lowest numbered!) the master of the region and
// increment the number of regions
regionMaster[nRegions] = own;
cellRegionMaster[nRegions] = own;
++nRegions;
}
else
@ -352,7 +356,7 @@ Foam::label Foam::removeFaces::compatibleRemoves
// See if owner becomes the master of the region (if its index
// is lower than the current master of the region)
regionMaster[region1] = min(own, regionMaster[region1]);
cellRegionMaster[region1] = min(own, cellRegionMaster[region1]);
}
}
else
@ -400,24 +404,29 @@ Foam::label Foam::removeFaces::compatibleRemoves
freedRegion = region0;
}
label master0 = regionMaster[region0];
label master1 = regionMaster[region1];
label master0 = cellRegionMaster[region0];
label master1 = cellRegionMaster[region1];
regionMaster[freedRegion] = -1;
regionMaster[keptRegion] = min(master0, master1);
cellRegionMaster[freedRegion] = -1;
cellRegionMaster[keptRegion] = min(master0, master1);
}
}
}
// Set size
regionMaster.setSize(nRegions);
cellRegionMaster.setSize(nRegions);
// Note:
// It is legal to have cellRegionMaster = -1 if the region has been created
// and then abandoned because it has been merged with another region
// HJ, 6/Sep/2019
// Various checks, additional scope for clarity and memory management
// - master is lowest numbered in any region
// - regions have more than 1 cell
{
// Number of cells per region
labelList nCells(regionMaster.size(), 0);
labelList nCells(cellRegionMaster.size(), 0);
// Loop through cell regions
forAll(cellRegion, cellI)
@ -430,15 +439,12 @@ Foam::label Foam::removeFaces::compatibleRemoves
// Region found, increment number of cells per region
++nCells[r];
if (cellI < regionMaster[r])
if (cellI < cellRegionMaster[r])
{
FatalErrorIn
(
"removeFaces::compatibleRemoves(const labelList&"
", labelList&, labelList&, labelList&)"
) << "Not lowest numbered! Cell: " << cellI
FatalErrorInFunction
<< "Not lowest numbered! Cell: " << cellI
<< " region: " << r
<< " region master: " << regionMaster[r]
<< " region master: " << cellRegionMaster[r]
<< abort(FatalError);
}
}
@ -449,11 +455,8 @@ Foam::label Foam::removeFaces::compatibleRemoves
{
if (nCells[regionI] == 1)
{
FatalErrorIn
(
"removeFaces::compatibleRemoves(const labelList&"
", labelList&, labelList&, labelList&)"
) << "Region " << regionI
FatalErrorInFunction
<< "Region " << regionI
<< " has only " << nCells[regionI] << " cell in it."
<< abort(FatalError);
}
@ -464,9 +467,9 @@ Foam::label Foam::removeFaces::compatibleRemoves
// Count number of used regions
label nUsedRegions = 0;
forAll(regionMaster, i)
forAll(cellRegionMaster, i)
{
if (regionMaster[i] != -1)
if (cellRegionMaster[i] != -1)
{
++nUsedRegions;
}

View file

@ -28,8 +28,6 @@ Description
Given list of faces to remove insert all the topology changes. Contains
helper function to get consistent set of faces to remove.
Not very well tested in parallel.
SourceFiles
removeFaces.C
@ -185,6 +183,7 @@ public:
// and always merge (if on same patch)
removeFaces(const polyMesh&, const scalar minCos);
// Member Functions
//- Given a set of faces to pierce calculates:
@ -210,7 +209,9 @@ public:
(
const labelList& piercedFaces,
const labelList& cellRegion,
const labelList& cellRegionMaster,
const labelList& pointRegionMaster,
labelList& cellRegionMaster,
TopoChangeEngine& ref
) const;

View file

@ -67,18 +67,9 @@ void Foam::removeFaces::mergeFaces
if (fp.edgeLoops().size() != 1)
{
writeOBJ(fp, mesh_.time().path()/"facesToBeMerged.obj");
FatalErrorIn
(
"template<class TopoChangeEngine>"
"void removeFaces::mergeFaces"
"\n("
"\n const labelList& cellRegion,"
"\n const labelList& cellRegionMaster,"
"\n const labelHashSet& pointsToRemove,"
"\n const labelList& faceLabels,"
"\n const TopoChangeEngine& ref"
"\n) const"
) << "Cannot merge faces " << faceLabels
FatalErrorInFunction
<< "Cannot merge faces " << faceLabels
<< " into single face since outside vertices " << fp.edgeLoops()
<< " do not form single loop but form " << fp.edgeLoops().size()
<< " loops instead." << abort(FatalError);
@ -131,18 +122,9 @@ void Foam::removeFaces::mergeFaces
if (masterIndex == -1)
{
writeOBJ(fp, mesh_.time().path()/"facesToBeMerged.obj");
FatalErrorIn
(
"template<class TopoChangeEngine>"
"void removeFaces::mergeFaces"
"\n("
"\n const labelList& cellRegion,"
"\n const labelList& cellRegionMaster,"
"\n const labelHashSet& pointsToRemove,"
"\n const labelList& faceLabels,"
"\n const TopoChangeEngine& ref"
"\n) const"
) << "Problem." << abort(FatalError);
FatalErrorInFunction
<< "Problem." << abort(FatalError);
}
@ -312,7 +294,8 @@ void Foam::removeFaces::setRefinement
(
const labelList& faceLabels,
const labelList& cellRegion,
const labelList& cellRegionMaster,
const labelList& pointRegionMaster,
labelList& cellRegionMaster,
TopoChangeEngine& ref
) const
@ -343,18 +326,8 @@ void Foam::removeFaces::setRefinement
if (!mesh_.isInternalFace(faceI))
{
FatalErrorIn
(
"template<class TopoChangeEngine>"
"\n"
"void removeFaces::setRefinement"
"\n("
"\n const labelList& faceLabels,"
"\n const labelList& cellRegion,"
"\n const labelList& cellRegionMaster,"
"\n TopoChangeEngine& ref"
"\n)"
) << "Face " << faceI << " is not an internal faces, therefore"
FatalErrorInFunction
<< "Face " << faceI << " is not an internal faces, therefore"
<< " it cannot be removed. Check faceLabels argument."
<< abort(FatalError);
}
@ -446,18 +419,8 @@ void Foam::removeFaces::setRefinement
mesh_.write();
// Write data for debugging
FatalErrorIn
(
"template<class TopoChangeEngine>"
"\n"
"void removeFaces::setRefinement"
"\n("
"\n const labelList& faceLabels,"
"\n const labelList& cellRegion,"
"\n const labelList& cellRegionMaster,"
"\n TopoChangeEngine& ref"
"\n)"
) << "Problem : edge has too few face neighbours:"
FatalErrorInFunction
<< "Problem : edge has too few face neighbours:"
<< eFaces << endl
<< "edge:" << edgeI
<< " vertices:" << e
@ -586,18 +549,8 @@ void Foam::removeFaces::setRefinement
const edge& e = mesh_.edges()[edgeI];
// Only found one boundary face. Problem
FatalErrorIn
(
"template<class TopoChangeEngine>"
"\n"
"void removeFaces::setRefinement"
"\n("
"\n const labelList& faceLabels,"
"\n const labelList& cellRegion,"
"\n const labelList& cellRegionMaster,"
"\n TopoChangeEngine& ref"
"\n)"
) << "Problem : edge would have one boundary face"
FatalErrorInFunction
<< "Problem : edge would have one boundary face"
<< " and one internal face using it." << endl
<< "Your remove pattern is probably incorrect." << endl
<< "edge:" << edgeI
@ -620,18 +573,7 @@ void Foam::removeFaces::setRefinement
{
const edge& e = mesh_.edges()[edgeI];
FatalErrorIn
(
"template<class TopoChangeEngine>"
"\n"
"void removeFaces::setRefinement"
"\n("
"\n const labelList& faceLabels,"
"\n const labelList& cellRegion,"
"\n const labelList& cellRegionMaster,"
"\n TopoChangeEngine& ref"
"\n)"
)
FatalErrorInFunction
<< "Problem : edge would get 1 face using it only"
<< " edge:" << edgeI
<< " nFaces:" << nFacesPerEdge[edgeI]
@ -779,18 +721,8 @@ void Foam::removeFaces::setRefinement
if (nRegion < 1)
{
FatalErrorIn
(
"template<class TopoChangeEngine>"
"\n"
"void removeFaces::setRefinement"
"\n("
"\n const labelList& faceLabels,"
"\n const labelList& cellRegion,"
"\n const labelList& cellRegionMaster,"
"\n TopoChangeEngine& ref"
"\n)"
) << "Problem with region number." << abort(FatalError);
FatalErrorInFunction
<< "Problem with region number." << abort(FatalError);
}
else if (nRegion == 1)
{
@ -837,18 +769,8 @@ void Foam::removeFaces::setRefinement
{
if (nbrRegion != myRegion)
{
FatalErrorIn
(
"template<class TopoChangeEngine>"
"\n"
"void removeFaces::setRefinement"
"\n("
"\n const labelList& faceLabels,"
"\n const labelList& cellRegion,"
"\n const labelList& cellRegionMaster,"
"\n TopoChangeEngine& ref"
"\n)"
) << "Inconsistent face region across coupled patches."
FatalErrorInFunction
<< "Inconsistent face region across coupled patches."
<< endl
<< "This side has for faceI:" << faceI
<< " region:" << myRegion << endl
@ -868,18 +790,8 @@ void Foam::removeFaces::setRefinement
// Second visit of this region
if (toNbrRegion[myRegion] != nbrRegion)
{
FatalErrorIn
(
"template<class TopoChangeEngine>"
"\n"
"void removeFaces::setRefinement"
"\n("
"\n const labelList& faceLabels,"
"\n const labelList& cellRegion,"
"\n const labelList& cellRegionMaster,"
"\n TopoChangeEngine& ref"
"\n)"
) << "Inconsistent face region across coupled patches."
FatalErrorInFunction
<< "Inconsistent face region across coupled patches."
<< endl
<< "This side has for faceI:" << faceI
<< " region:" << myRegion
@ -938,18 +850,8 @@ void Foam::removeFaces::setRefinement
{
if (nEdgesPerPoint[pointI] == 1)
{
FatalErrorIn
(
"template<class TopoChangeEngine>"
"\n"
"void removeFaces::setRefinement"
"\n("
"\n const labelList& faceLabels,"
"\n const labelList& cellRegion,"
"\n const labelList& cellRegionMaster,"
"\n TopoChangeEngine& ref"
"\n)"
) << "Problem : point would get 1 edge using it only."
FatalErrorInFunction
<< "Problem : point would get 1 edge using it only."
<< " pointI:" << pointI
<< " coord:" << mesh_.points()[pointI]
<< abort(FatalError);
@ -1060,11 +962,43 @@ void Foam::removeFaces::setRefinement
ref.setAction(polyRemovePoint(pointI, -1));
}
// Add master cells for correct mapping
forAll (cellRegionMaster, regionI)
{
// Note: it is legal to have cellRegionMaster = -1 if the region
// has been created and them marged into another region.
// Such masters will also have pointRegionMaster = -1 and should
// be ignored. HJ, 6/Sep/2019
// Additionally protect for old directTopoChangers which do not
// identify points for mapping. Non-existent pointRegionMaster
// is rejected
if (cellRegionMaster[regionI] > -1 && pointRegionMaster[regionI] > -1)
{
// Add master cell from master point for correct mapping
cellRegionMaster[regionI] =
ref.setAction
(
polyAddCell
(
pointRegionMaster[regionI], // masterPointID
-1, // masterEdgeID
-1, // masterFaceID
-1, // masterCellID
mesh_.cellZones().whichZone(cellRegionMaster[regionI])
)
);
}
}
// Remove cells
forAll(cellRegion, cellI)
{
label region = cellRegion[cellI];
// Old check is acceptable: for mapping from point, the cellRegionMaster
// has been replaced in polyAddPoint
// HJ, 6/Sep/2019
if (region != -1 && (cellI != cellRegionMaster[region]))
{
ref.setAction(polyRemoveCell(cellI, cellRegionMaster[region]));
@ -1087,18 +1021,8 @@ void Foam::removeFaces::setRefinement
if (rFaces.size() <= 1)
{
FatalErrorIn
(
"template<class TopoChangeEngine>"
"\n"
"void removeFaces::setRefinement"
"\n("
"\n const labelList& faceLabels,"
"\n const labelList& cellRegion,"
"\n const labelList& cellRegionMaster,"
"\n TopoChangeEngine& ref"
"\n)"
) << "Region: " << regionI
FatalErrorInFunction
<< "Region: " << regionI
<< " contains only these faces: " << rFaces
<< abort(FatalError);
}

View file

@ -634,10 +634,16 @@ Foam::labelList Foam::undoableMeshCutter::removeSplitFaces
// Insert all commands to combine cells. Never fails so don't have to
// test for success.
// Added dummy pointRegionMaster, which will be ignored.
// HJ, 6/Sep/2019
labelList pointRegionMaster(cellRegionMaster.size(), label(-1));
faceRemover().setRefinement
(
facesToRemove,
cellRegion,
pointRegionMaster,
cellRegionMaster,
meshMod
);

View file

@ -1260,10 +1260,8 @@ void Foam::polyhedralRefinement::setRefinementInstruction
// cell points
const labelList& cPoints = meshCellPoints[cellI];
FatalErrorIn
(
"polyhedralRefinement::setRefinementInstruction(...)"
) << "Cell " << cellI
FatalErrorInFunction
<< "Cell " << cellI
<< " of level " << cellLevel_[cellI]
<< " does not seem to have enough points of "
<< " lower level" << endl
@ -1967,11 +1965,73 @@ void Foam::polyhedralRefinement::setUnrefinementInstruction
}
}
// Find point region master for every cell region. This is the central point
// from which the coarse cell will be made
// The property of the point region master is that all cells that touch it
// have the same cell region index
// HJ, 6/Sep/2019
labelList pointRegionMaster(cellRegionMaster.size(), label(-1));
// Get point-cell addressing
const labelListList& pc = mesh_.pointCells();
forAll (splitPointsToUnrefine_, i)
{
const labelList& curPc = pc[splitPointsToUnrefine_[i]];
label curRegion = -1;
forAll (curPc, curPcI)
{
if (curRegion == -1)
{
// First region found. Grab it
curRegion = cellRegion[curPc[curPcI]];
}
else
{
// Region already found. Check that all other cells that
// touch this point have the same region
if (curRegion != cellRegion[curPc[curPcI]])
{
// Error: different region cells touching in split point
// This is not a valid unrefinement pattern
FatalErrorIn
(
"polyhedralRefinement::setUnrefinementInstruction"
"(polyTopoChange& ref)"
) << "Different region cells touching in split point."
<< abort(FatalError);
}
}
}
// Record point region master
if (curRegion > -1)
{
pointRegionMaster[curRegion] = splitPointsToUnrefine_[i];
}
else
{
// Error: Cannot find region for point
FatalErrorIn
(
"polyhedralRefinement::setUnrefinementInstruction"
"(polyTopoChange& ref)"
) << "Different region cells touching in split point."
<< abort(FatalError);
}
}
Info<< "pointRegionMaster: " << pointRegionMaster << nl
<< " cellRegionMaster: " << cellRegionMaster << endl;
// Insert all commands to combine cells
faceRemover_.setRefinement
(
facesToRemove,
cellRegion,
pointRegionMaster,
cellRegionMaster,
ref
);

View file

@ -2329,7 +2329,7 @@ void Foam::prismatic2DRefinement::setRefinementInstruction
if (minPointI != labelMax && minPointI != mesh_.nPoints())
{
FatalErrorIn("prismatic2DRefinement::setRefinementInstruction(...)")
FatalErrorInFunction
<< "Added point labels not consecutive to existing mesh points."
<< nl
<< "mesh_.nPoints():" << mesh_.nPoints()
@ -2352,10 +2352,8 @@ void Foam::prismatic2DRefinement::setUnrefinementInstruction
// Check whether the refinementLevelIndicator is valid
if (refinementLevelIndicator_.size() != mesh_.nCells())
{
FatalErrorIn
(
"prismatic2DRefinement::setUnrefinementInstruction(...)"
) << "Refinement level indicator list has invalid size: "
FatalErrorInFunction
<< "Refinement level indicator list has invalid size: "
<< refinementLevelIndicator_.size()
<< ", number of cells: " << mesh_.nCells()
<< nl
@ -2379,11 +2377,8 @@ void Foam::prismatic2DRefinement::setUnrefinementInstruction
{
if (cellLevel_[cellI] < 0)
{
FatalErrorIn
(
"prismatic2DRefinement::setUnrefinementInstruction"
"(polyTopoChange& ref)"
) << "Illegal cell level " << cellLevel_[cellI]
FatalErrorInFunction
<< "Illegal cell level " << cellLevel_[cellI]
<< " for cell " << cellI
<< abort(FatalError);
}
@ -2526,11 +2521,70 @@ void Foam::prismatic2DRefinement::setUnrefinementInstruction
}
}
// Find point region master for every cell region. This is the central point
// from which the coarse cell will be made
// The property of the point region master is that all cells that touch it
// have the same cell region index
// HJ, 6/Sep/2019
labelList pointRegionMaster(cellRegionMaster.size(), label(-1));
// Get point-cell addressing
const labelListList& pc = mesh_.pointCells();
forAll (splitPointsToUnrefine_, i)
{
const labelList& curPc = pc[splitPointsToUnrefine_[i]];
label curRegion = -1;
forAll (curPc, curPcI)
{
if (curRegion == -1)
{
// First region found. Grab it
curRegion = cellRegion[curPc[curPcI]];
}
else
{
// Region already found. Check that all other cells that
// touch this point have the same region
if (curRegion != cellRegion[curPc[curPcI]])
{
// Error: different region cells touching in split point
// This is not a valid unrefinement pattern
FatalErrorIn
(
"polyhedralRefinement::setUnrefinementInstruction"
"(polyTopoChange& ref)"
) << "Different region cells touching in split point."
<< abort(FatalError);
}
}
}
// Record point region master
if (curRegion > -1)
{
pointRegionMaster[curRegion] = splitPointsToUnrefine_[i];
}
else
{
// Error: Cannot find region for point
FatalErrorIn
(
"polyhedralRefinement::setUnrefinementInstruction"
"(polyTopoChange& ref)"
) << "Different region cells touching in split point."
<< abort(FatalError);
}
}
// Insert all commands to combine cells
faceRemover_.setRefinement
(
facesToRemove,
cellRegion,
pointRegionMaster,
cellRegionMaster,
ref
);

View file

@ -23,6 +23,7 @@ License
Author
Vuko Vukcevic, Wikki Ltd. All rights reserved.
Hrvoje Jasak, Wikki Ltd.
Notes
Generalisation of hexRef8 for polyhedral cells and refactorisation into mesh
@ -423,6 +424,54 @@ Foam::label Foam::refinement::countAnchors
return nAnchors;
}
void Foam::refinement::adjustRefLevel
(
label& curNewCellLevel,
const label oldCellI
)
{
if (oldCellI == -1)
{
// This cell is inflated (does not originate from other cell), set
// cell level to -1
curNewCellLevel = -1;
}
else
{
// This cell has either been added based on another cell or it
// hasn't changed. Update new cell level according to refinement
// level indicator and old cell level
// Get refinement status of the old cell
const label& refStatus = refinementLevelIndicator_[oldCellI];
if (refStatus == UNREFINED)
{
// New cell has been obtained by unrefining other cells - this
// is the remaining "master" cell. Decrement cell level
curNewCellLevel = cellLevel_[oldCellI] - 1;
}
else if (refStatus == UNCHANGED)
{
// Cell hasn't been changed during this refinement, copy old
// cell level
curNewCellLevel = cellLevel_[oldCellI];
}
else if (refStatus == REFINED)
{
// Cell has been refined, increment cell level
curNewCellLevel = cellLevel_[oldCellI] + 1;
}
else
{
FatalErrorInFunction
<< "Invalid refinement status detected: "
<< refStatus << nl
<< "Old cell index: " << oldCellI << abort(FatalError);
}
}
}
void Foam::refinement::checkInternalOrientation
(
@ -1386,51 +1435,21 @@ void Foam::refinement::updateMesh(const mapPolyMesh& map)
// Loop through all new cells
forAll (cellMap, newCellI)
{
// Get index of the corresponding old cell
const label& oldCellI = cellMap[newCellI];
adjustRefLevel(newCellLevel[newCellI], cellMap[newCellI]);
}
if (oldCellI == -1)
{
// This cell is inflated (does not originate from other cell), set
// cell level to -1
newCellLevel[newCellI] = -1;
}
else
{
// This cell has either been added based on another cell or it
// hasn't changed. Update new cell level according to refinement
// level indicator and old cell level
// Do cells from points: unrefinement
// Note: should other types be done as well?
// HJ, 9/Sep/2019
const List<objectMap>& cellsFromPoints = map.cellsFromPointsMap();
// Get refinement status of the old cell
const label& refStatus = refinementLevelIndicator_[oldCellI];
if (refStatus == UNREFINED)
{
// New cell has been obtained by unrefining other cells - this
// is the remaining "master" cell. Decrement cell level
newCellLevel[newCellI] = cellLevel_[oldCellI] - 1;
}
else if (refStatus == UNCHANGED)
{
// Cell hasn't been changed during this refinement, copy old
// cell level
newCellLevel[newCellI] = cellLevel_[oldCellI];
}
else if (refStatus == REFINED)
{
// Cell has been refined, increment cell level
newCellLevel[newCellI] = cellLevel_[oldCellI] + 1;
}
else
{
FatalErrorIn
(
"refinement::updateMesh(const mapPolyMesh& map)"
) << "Invalid refinement status detected: " << refStatus << nl
<< "Old cell index: " << oldCellI << nl
<< "New cell index: " << newCellI << abort(FatalError);
}
}
forAll (cellsFromPoints, cfpI)
{
adjustRefLevel
(
newCellLevel[cellsFromPoints[cfpI].index()],
cellsFromPoints[cfpI].masterObjects()[0]
);
}
// Transfer the new cell level into the data member

View file

@ -232,6 +232,13 @@ protected:
const label anchorLevel
) const;
//- Adjust cell refinement level after topo change
void adjustRefLevel
(
label& curNewCellLevel,
const label oldCellI
);
// Debug functions