Commit 774db288 authored by mattijs's avatar mattijs
Browse files

ENH: regionToCell: allow multiple regions and erosion

parent 100a25ee
......@@ -167,8 +167,10 @@ FoamFile
// source regionToCell;
// sourceInfo
// {
// set c0; // name of cellSet giving mesh subset
// insidePoint (1 2 3); // point inside region to select
// set c0; // optional name of cellSet giving mesh subset
// insidePoints ((1 2 3)); // points inside region to select
// nErode 0; // optional number of layers to erode
// // selection
// }
//
// // Cells underneath plane such that volume is reached. E.g. for use
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -24,12 +24,10 @@ License
\*---------------------------------------------------------------------------*/
#include "regionToCell.H"
#include "polyMesh.H"
#include "regionSplit.H"
#include "globalMeshData.H"
#include "emptyPolyPatch.H"
#include "cellSet.H"
#include "syncTools.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
......@@ -49,104 +47,339 @@ addToRunTimeSelectionTable(topoSetSource, regionToCell, istream);
Foam::topoSetSource::addToUsageTable Foam::regionToCell::usage_
(
regionToCell::typeName,
"\n Usage: regionToCell subCellSet (x y z)\n\n"
" Select all cells in the connected region containing point.\n"
" If started inside the subCellSet keeps to it;\n"
" if started outside stays outside.\n"
"\n Usage: regionToCell subCellSet (pt0 .. ptn)\n\n"
" Select all cells in the connected region containing"
" points (pt0..ptn).\n"
);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::regionToCell::combine(topoSet& set, const bool add) const
void Foam::regionToCell::markRegionFaces
(
const boolList& selectedCell,
boolList& regionFace
) const
{
label cellI = mesh_.findCell(insidePoint_);
// Internal faces
const labelList& faceOwner = mesh_.faceOwner();
const labelList& faceNeighbour = mesh_.faceNeighbour();
forAll(faceNeighbour, faceI)
{
if
(
selectedCell[faceOwner[faceI]]
!= selectedCell[faceNeighbour[faceI]]
)
{
regionFace[faceI] = true;
}
}
// Swap neighbour selectedCell state
boolList nbrSelected;
syncTools::swapBoundaryCellList(mesh_, selectedCell, nbrSelected);
// Boundary faces
const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
forAll(pbm, patchI)
{
const polyPatch& pp = pbm[patchI];
const labelUList& faceCells = pp.faceCells();
forAll(faceCells, i)
{
label faceI = pp.start()+i;
label bFaceI = faceI-mesh_.nInternalFaces();
if
(
selectedCell[faceCells[i]]
!= selectedCell[nbrSelected[bFaceI]]
)
{
regionFace[faceI] = true;
}
}
}
}
Foam::boolList Foam::regionToCell::findRegions
(
const bool verbose,
const regionSplit& cellRegion
) const
{
boolList keepRegion(cellRegion.nRegions(), false);
forAll(insidePoints_, i)
{
// Find the region containing the insidePoint
label cellI = mesh_.findCell(insidePoints_[i]);
label keepRegionI = -1;
label keepProcI = -1;
if (cellI != -1)
{
keepRegionI = cellRegion[cellI];
keepProcI = Pstream::myProcNo();
}
reduce(keepRegionI, maxOp<label>());
keepRegion[keepRegionI] = true;
reduce(keepProcI, maxOp<label>());
if (keepProcI == -1)
{
FatalErrorIn
(
"outsideCellSelection::findRegions"
"(const bool, const regionSplit&)"
) << "Did not find " << insidePoints_[i]
<< " in mesh." << " Mesh bounds are " << mesh_.bounds()
<< exit(FatalError);
}
if (verbose)
{
Info<< "Found location " << insidePoints_[i]
<< " in cell " << cellI << " on processor " << keepProcI
<< " in global region " << keepRegionI
<< " out of " << cellRegion.nRegions() << " regions." << endl;
}
}
return keepRegion;
}
// Load the subset of cells
void Foam::regionToCell::unselectOutsideRegions
(
boolList& selectedCell
) const
{
// Determine faces on the edge of selectedCell
boolList blockedFace(mesh_.nFaces(), false);
markRegionFaces(selectedCell, blockedFace);
// Determine regions
regionSplit cellRegion(mesh_, blockedFace);
// Determine regions containing insidePoints_
boolList keepRegion(findRegions(true, cellRegion));
// Go back to bool per cell
forAll(cellRegion, cellI)
{
Info<< "Loading subset " << setName_ << " to delimit search region."
<< endl;
cellSet subSet(mesh_, setName_);
if (!keepRegion[cellRegion[cellI]])
{
selectedCell[cellI] = false;
}
}
}
boolList inSubset(mesh_.nCells(), false);
forAllConstIter(cellSet, subSet, iter)
void Foam::regionToCell::shrinkRegions
(
boolList& selectedCell
) const
{
// Select points on unselected cells and boundary
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
boolList boundaryPoint(mesh_.nPoints(), false);
const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
forAll(pbm, patchI)
{
const polyPatch& pp = pbm[patchI];
if (!pp.coupled() && !isA<emptyPolyPatch>(pp))
{
inSubset[iter.key()] = true;
forAll(pp, i)
{
const face& f = pp[i];
forAll(f, fp)
{
boundaryPoint[f[fp]] = true;
}
}
}
}
if (cellI != -1 && inSubset[cellI])
forAll(selectedCell, cellI)
{
if (!selectedCell[cellI])
{
Pout<< "Point " << insidePoint_ << " is inside cellSet "
<< setName_ << endl
<< "Collecting all cells connected to " << cellI
<< " and inside cellSet " << setName_ << endl;
const labelList& cPoints = mesh_.cellPoints(cellI);
forAll(cPoints, i)
{
boundaryPoint[cPoints[i]] = true;
}
}
else
}
syncTools::syncPointList(mesh_, boundaryPoint, orEqOp<bool>(), false);
// Select all cells using these points
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
label nChanged = 0;
forAll(boundaryPoint, pointI)
{
if (boundaryPoint[pointI])
{
Pout<< "Point " << insidePoint_ << " is outside cellSet "
<< setName_ << endl
<< "Collecting all cells connected to " << cellI
<< " and outside cellSet " << setName_ << endl;
const labelList& pCells = mesh_.pointCells(pointI);
forAll(pCells, i)
{
label cellI = pCells[i];
if (selectedCell[cellI])
{
selectedCell[cellI] = false;
nChanged++;
}
}
}
}
}
void Foam::regionToCell::erode
(
boolList& selectedCell
) const
{
//Info<< "Entering shrinkRegions:" << count(selectedCell) << endl;
//generateField("selectedCell_before", selectedCell)().write();
// Now erode and see which regions get disconnected
boolList shrunkSelectedCell(selectedCell);
for (label iter = 0; iter < nErode_; iter++)
{
shrinkRegions(shrunkSelectedCell);
}
//Info<< "After shrinking:" << count(shrunkSelectedCell) << endl;
//generateField("shrunkSelectedCell", shrunkSelectedCell)().write();
// Determine faces on the edge of shrunkSelectedCell
boolList blockedFace(mesh_.nFaces(), false);
markRegionFaces(shrunkSelectedCell, blockedFace);
// Find disconnected regions
regionSplit cellRegion(mesh_, blockedFace);
// Determine regions containing insidePoints
boolList keepRegion(findRegions(true, cellRegion));
// Get coupled cell status
label nInt = mesh_.nInternalFaces();
boolList neiSet(mesh_.nFaces()-nInt, false);
for (label faceI = nInt; faceI < mesh_.nFaces(); faceI++)
// Extract cells in regions that are not to be kept.
boolList removeCell(mesh_.nCells(), false);
forAll(cellRegion, cellI)
{
if (shrunkSelectedCell[cellI] && !keepRegion[cellRegion[cellI]])
{
neiSet[faceI-nInt] = inSubset[mesh_.faceOwner()[faceI]];
removeCell[cellI] = true;
}
syncTools::swapBoundaryFaceList(mesh_, neiSet);
}
//Info<< "removeCell before:" << count(removeCell) << endl;
//generateField("removeCell_before", removeCell)().write();
// Find faces inbetween subSet and non-subset.
for (label faceI = 0; faceI < nInt; faceI++)
// Grow removeCell
for (label iter = 0; iter < nErode_; iter++)
{
// Grow selected cell in regions that are not for keeping
boolList boundaryPoint(mesh_.nPoints(), false);
forAll(removeCell, cellI)
{
bool ownInSet = inSubset[mesh_.faceOwner()[faceI]];
bool neiInSet = inSubset[mesh_.faceNeighbour()[faceI]];
blockedFace[faceI] = (ownInSet != neiInSet);
if (removeCell[cellI])
{
const labelList& cPoints = mesh_.cellPoints(cellI);
forAll(cPoints, i)
{
boundaryPoint[cPoints[i]] = true;
}
}
}
for (label faceI = nInt; faceI < mesh_.nFaces(); faceI++)
syncTools::syncPointList(mesh_, boundaryPoint, orEqOp<bool>(), false);
// Select all cells using these points
label nChanged = 0;
forAll(boundaryPoint, pointI)
{
bool ownInSet = inSubset[mesh_.faceOwner()[faceI]];
bool neiInSet = neiSet[faceI-nInt];
blockedFace[faceI] = (ownInSet != neiInSet);
if (boundaryPoint[pointI])
{
const labelList& pCells = mesh_.pointCells(pointI);
forAll(pCells, i)
{
label cellI = pCells[i];
if (!removeCell[cellI])
{
removeCell[cellI] = true;
nChanged++;
}
}
}
}
}
// Find connected regions without crossing boundary of the cellset.
regionSplit regions(mesh_, blockedFace);
//Info<< "removeCell after:" << count(removeCell) << endl;
//generateField("removeCell_after", removeCell)().write();
// Get the region containing the insidePoint
label regionI = -1;
if (cellI != -1)
// Unmark removeCell
forAll(removeCell, cellI)
{
// On processor that has found cell.
regionI = regions[cellI];
if (removeCell[cellI])
{
selectedCell[cellI] = false;
}
}
}
reduce(regionI, maxOp<label>());
void Foam::regionToCell::combine(topoSet& set, const bool add) const
{
// Note: wip. Select cells first
boolList selectedCell(mesh_.nCells(), true);
if (regionI == -1)
if (setName_.size() && setName_ != "none")
{
WarningIn
(
"regionToCell::combine(topoSet&, const bool) const"
) << "Point " << insidePoint_
<< " is not inside the mesh." << nl
<< "Bounding box of the mesh:" << mesh_.bounds()
Info<< "Loading subset " << setName_ << " to delimit search region."
<< endl;
return;
cellSet subSet(mesh_, setName_);
selectedCell = false;
forAllConstIter(cellSet, subSet, iter)
{
selectedCell[iter.key()] = true;
}
}
// Pick up the cells of the region
const labelList regionCells(findIndices(regions, regionI));
unselectOutsideRegions(selectedCell);
forAll(regionCells, i)
if (nErode_ > 0)
{
addOrDelete(set, regionCells[i], add);
erode(selectedCell);
}
forAll(selectedCell, cellI)
{
if (selectedCell[cellI])
{
addOrDelete(set, cellI, add);
}
}
}
......@@ -158,12 +391,14 @@ Foam::regionToCell::regionToCell
(
const polyMesh& mesh,
const word& setName,
const point& insidePoint
const pointField& insidePoints,
const label nErode
)
:
topoSetSource(mesh),
setName_(setName),
insidePoint_(insidePoint)
insidePoints_(insidePoints),
nErode_(nErode)
{}
......@@ -175,8 +410,14 @@ Foam::regionToCell::regionToCell
)
:
topoSetSource(mesh),
setName_(dict.lookup("set")),
insidePoint_(dict.lookup("insidePoint"))
setName_(dict.lookupOrDefault<word>("set", "none")),
insidePoints_
(
dict.found("insidePoints")
? dict.lookup("insidePoints")
: dict.lookup("insidePoint")
),
nErode_(dict.lookupOrDefault("nErode", 0))
{}
......@@ -189,7 +430,8 @@ Foam::regionToCell::regionToCell
:
topoSetSource(mesh),
setName_(checkIs(is)),
insidePoint_(checkIs(is))
insidePoints_(checkIs(is)),
nErode_(0)
{}
......@@ -209,15 +451,15 @@ void Foam::regionToCell::applyToSet
{
if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
{
Info<< " Adding all cells of connected region containing point "
<< insidePoint_ << " ..." << endl;
Info<< " Adding all cells of connected region containing points "
<< insidePoints_ << " ..." << endl;
combine(set, true);
}
else if (action == topoSetSource::DELETE)
{
Info<< " Removing all cells of connected region containing point "
<< insidePoint_ << " ..." << endl;
Info<< " Removing all cells of connected region containing points "
<< insidePoints_ << " ..." << endl;
combine(set, false);
}
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -26,7 +26,20 @@ Class
Description
TopoSetSource. Select cells belonging to topological connected region
(that contains given point)
(that contains given points)
In dictionary input:
// optional name of cellSet delimiting search
set c0;
//- Number of cell layers to erode mesh to detect holes in the mesh
// Set to 0 if not used.
nErode 3;
// points inside region to select
insidePoints ((1 2 3));
SourceFiles
regionToCell.C
......@@ -37,12 +50,15 @@ SourceFiles
#define regionToCell_H
#include "topoSetSource.H"
#include "boolList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class regionSplit;
/*---------------------------------------------------------------------------*\
Class regionToCell Declaration
\*---------------------------------------------------------------------------*/
......@@ -58,14 +74,38 @@ class regionToCell
static addToUsageTable usage_;
//- Name of cellSet to keep to
word setName_;
const word setName_;
//- Coordinate(s) that is inside connected region
const pointField insidePoints_;
//- Coordinate that is inside connected region
point insidePoint_;
//- Number of layers to erode
const label nErode_;
// Private Member Functions
//- Mark faces inbetween selected and unselected elements
void markRegionFaces
(
const boolList& selectedCell,
boolList& regionFace
) const;
//- Determine for every disconnected region in the mesh whether
// it contains a locationInMesh
boolList findRegions(const bool verbose, const regionSplit&) const;
//- Unselect regions not containing a locationInMesh
void unselectOutsideRegions(boolList& selectedCell) const;
//- Unselect one layer of cells from selectedCell
void shrinkRegions(boolList& selectedCell) const;
//- Erode a given number of layers from selectedCell. Remove any
// region that gets disconnected that way.
void erode(boolList& selectedCell) const;
void combine(topoSet& set, const bool add) const;
......@@ -81,7 +121,8 @@ public:
(
const polyMesh& mesh,
const word& setName,
const point& insidePoint
const pointField& insidePoints,
const label nErode
);
//- Construct from dictionary
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment