Commit 9973c378 authored by Mark OLESEN's avatar Mark OLESEN
Browse files

ENH: refactor cutting-plane cell selection

- avoid duplicate code by relocating cellZone selection and bounding box
  sub-selection into cuttingPlane and cuttingSurfaceBaseSelection.
  Allows reuse by inherited classes (sampledPlane, surfMeshSamplePlane).
parent a8ef9e97
......@@ -21,7 +21,9 @@ sampledSet/shortestPath/shortestPathSet.C
surface/cutting/cuttingPlane.C
surface/cutting/cuttingPlaneCuts.C
surface/cutting/cuttingPlaneSelection.C
surface/cutting/cuttingSurfaceBase.C
surface/cutting/cuttingSurfaceBaseSelection.C
surface/distanceSurface/distanceSurface.C
surface/isoSurface/isoSurface.C
surface/isoSurface/isoSurfaceCell.C
......
......@@ -47,153 +47,16 @@ namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::sampledPlane::checkBoundsIntersection
(
const plane& pln,
const boundBox& meshBb
) const
{
// Verify specified bounding box
if (!bounds_.empty())
{
// Bounding box does not overlap with (global) mesh!
if (!bounds_.overlaps(meshBb))
{
WarningInFunction
<< nl
<< name() << " : "
<< "Bounds " << bounds_
<< " do not overlap the mesh bounding box " << meshBb
<< nl << endl;
}
// Plane does not intersect the bounding box
if (!bounds_.intersects(pln))
{
WarningInFunction
<< nl
<< name() << " : "
<< "Plane "<< pln << " does not intersect the bounds "
<< bounds_
<< nl << endl;
}
}
// Plane does not intersect the (global) mesh!
if (!meshBb.intersects(pln))
{
WarningInFunction
<< nl
<< name() << " : "
<< "Plane "<< pln << " does not intersect the mesh bounds "
<< meshBb
<< nl << endl;
}
}
Foam::bitSet Foam::sampledPlane::cellSelection
(
const bool warnIntersect
) const
Foam::bitSet Foam::sampledPlane::cellSelection(const bool warn) const
{
// Zones requested and in use?
const bool hasZones =
returnReduce
(
(-1 != mesh().cellZones().findIndex(zoneNames_)),
andOp<bool>()
);
bitSet cellsToSelect;
if (hasZones)
{
cellsToSelect = mesh().cellZones().selection(zoneNames_);
}
// Subset the zoned cells with the bounds_.
// For a full mesh, use the bounds to define the cell selection.
// If there are zones cells, use them to build the effective mesh
// bound box.
// Note that for convenience we use cell centres here instead of the
// cell points, since it will only be used for checking.
boundBox meshBb;
if (bounds_.empty())
{
// No bounds restriction, but will need the effective mesh boundBox
// for checking intersections
if (hasZones && warnIntersect)
{
const auto& cellCentres = static_cast<const fvMesh&>(mesh()).C();
for (const label celli : cellsToSelect)
{
const point& cc = cellCentres[celli];
meshBb.add(cc);
}
meshBb.reduce();
}
else
{
meshBb = mesh().bounds(); // use the regular mesh bounding box
}
}
else
{
const auto& cellCentres = static_cast<const fvMesh&>(mesh()).C();
// Only examine cells already set
if (hasZones)
{
for (const label celli : cellsToSelect)
{
const point& cc = cellCentres[celli];
meshBb.add(cc);
if (!bounds_.contains(cc))
{
cellsToSelect.unset(celli);
}
}
meshBb.reduce();
}
else
{
const label len = mesh().nCells();
cellsToSelect.resize(len);
for (label celli=0; celli < len; ++celli)
{
const point& cc = cellCentres[celli];
if (bounds_.contains(cc))
{
cellsToSelect.set(celli);
}
}
meshBb = mesh().bounds(); // use the regular mesh bounding box
}
}
if (warnIntersect)
{
checkBoundsIntersection(*this, meshBb);
}
return cellsToSelect;
return cuttingPlane::cellSelection
(
mesh(),
bounds_,
zoneNames_,
name(),
warn
);
}
......@@ -257,9 +120,10 @@ Foam::sampledPlane::sampledPlane
if (dict.found("coordinateSystem"))
{
coordinateSystem cs(mesh, dict.subDict("coordinateSystem"));
plane& pln = planeDesc();
const point orig = cs.globalPosition(planeDesc().origin());
const vector norm = cs.globalVector(planeDesc().normal());
const point orig = cs.globalPosition(pln.origin());
const vector norm = cs.globalVector(pln.normal());
if (debug)
{
......@@ -269,8 +133,8 @@ Foam::sampledPlane::sampledPlane
<< " defined within a local coordinateSystem" << endl;
}
// Assign the plane description
static_cast<plane&>(*this) = plane(orig, norm);
// Reassign the plane
pln = plane(orig, norm);
}
......@@ -331,7 +195,7 @@ bool Foam::sampledPlane::update()
sampledSurface::clearGeom();
performCut(mesh(), triangulate_, this->cellSelection(true));
performCut(mesh(), triangulate_, cellSelection(true));
if (debug)
{
......
......@@ -56,7 +56,7 @@ Usage
bounds | limit with bounding box | no |
zone | limit to cell zone (name or regex) | no |
zones | limit to cell zones (names, regexs) | no |
coordinateSystem | plane relative to given coordinate system | no |
coordinateSystem | define plane within given coordinate system | no |
\endtable
SeeAlso
......@@ -98,10 +98,10 @@ class sampledPlane
//- The zone or zones in which cutting is to occur
wordRes zoneNames_;
//- Optional bounding box to trim triangles against
//- Optional bounding box to trim against
const boundBox bounds_;
//- Triangulated faces or keep faces as is
//- Triangulate faces or not
const bool triangulate_;
//- Track if the surface needs an update
......@@ -110,18 +110,11 @@ class sampledPlane
// Private Member Functions
//- Check and warn if bounding box does not intersect mesh or plane
void checkBoundsIntersection
(
const plane& pln,
const boundBox& meshBb
) const;
//- Define cell selection from zones and bounding box.
// Optionally check and warn if the plane does not intersect
// with the bounds of the mesh (or submesh) or if the bounding box
// does not overlap with the mesh (or submesh)
bitSet cellSelection(const bool warnIntersect=false) const;
bitSet cellSelection(const bool warn=false) const;
//- Sample volume field onto surface faces
......
......@@ -47,152 +47,16 @@ namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::surfMeshSamplePlane::checkBoundsIntersection
(
const plane& pln,
const boundBox& meshBb
) const
{
// Verify specified bounding box
if (!bounds_.empty())
{
// Bounding box does not overlap with (global) mesh!
if (!bounds_.overlaps(meshBb))
{
WarningInFunction
<< nl
<< name() << " : "
<< "Bounds " << bounds_
<< " do not overlap the mesh bounding box " << mesh().bounds()
<< nl << endl;
}
// Plane does not intersect the bounding box
if (!bounds_.intersects(pln))
{
WarningInFunction
<< nl
<< name() << " : "
<< "Plane "<< pln << " does not intersect the bounds "
<< bounds_
<< nl << endl;
}
}
// Plane does not intersect the (global) mesh!
if (!meshBb.intersects(pln))
{
WarningInFunction
<< nl
<< name() << " : "
<< "Plane "<< pln << " does not intersect the mesh bounds "
<< meshBb
<< nl << endl;
}
}
Foam::bitSet Foam::surfMeshSamplePlane::cellSelection
(
const bool warnIntersect
) const
Foam::bitSet Foam::surfMeshSamplePlane::cellSelection(const bool warn) const
{
// Zones requested and in use?
const bool hasZones =
returnReduce
(
(-1 != mesh().cellZones().findIndex(zoneNames_)),
andOp<bool>()
);
bitSet cellsToSelect;
if (hasZones)
{
cellsToSelect = mesh().cellZones().selection(zoneNames_);
}
// Subset the zoned cells with the bounds_.
// For a full mesh, use the bounds to define the cell selection.
// If there are zones cells, use them to build the effective mesh
// bound box.
// Note that for convenience we use cell centres here instead of the
// cell points, since it will only be used for checking.
boundBox meshBb;
if (bounds_.empty())
{
// No bounds restriction, but will need the effective mesh boundBox
// for checking intersections
if (hasZones && warnIntersect)
{
const auto& cellCentres = static_cast<const fvMesh&>(mesh()).C();
for (const label celli : cellsToSelect)
{
const point& cc = cellCentres[celli];
meshBb.add(cc);
}
meshBb.reduce();
}
else
{
meshBb = mesh().bounds(); // use the regular mesh bounding box
}
}
else
{
const auto& cellCentres = static_cast<const fvMesh&>(mesh()).C();
// Only examine cells already set
if (hasZones)
{
for (const label celli : cellsToSelect)
{
const point& cc = cellCentres[celli];
meshBb.add(cc);
if (!bounds_.contains(cc))
{
cellsToSelect.unset(celli);
}
}
meshBb.reduce();
}
else
{
const label len = mesh().nCells();
cellsToSelect.resize(len);
for (label celli=0; celli < len; ++celli)
{
const point& cc = cellCentres[celli];
if (bounds_.contains(cc))
{
cellsToSelect.set(celli);
}
}
meshBb = mesh().bounds(); // use the regular mesh bounding box
}
}
if (warnIntersect)
{
checkBoundsIntersection(*this, meshBb);
}
return cellsToSelect;
return cuttingPlane::cellSelection
(
mesh(),
bounds_,
zoneNames_,
name(),
warn
);
}
......@@ -256,9 +120,10 @@ Foam::surfMeshSamplePlane::surfMeshSamplePlane
if (dict.found("coordinateSystem"))
{
coordinateSystem cs(mesh, dict.subDict("coordinateSystem"));
plane& pln = planeDesc();
const point orig = cs.globalPosition(planeDesc().origin());
const vector norm = cs.globalVector(planeDesc().normal());
const point orig = cs.globalPosition(pln.origin());
const vector norm = cs.globalVector(pln.normal());
if (debug)
{
......@@ -268,8 +133,8 @@ Foam::surfMeshSamplePlane::surfMeshSamplePlane
<< " defined within a local coordinateSystem" << endl;
}
// Assign the plane description
static_cast<plane&>(*this) = plane(orig, norm);
// Reassign the plane
pln = plane(orig, norm);
}
if (debug)
......@@ -325,7 +190,7 @@ bool Foam::surfMeshSamplePlane::update()
return false;
}
performCut(mesh(), triangulate_, this->cellSelection(true));
performCut(mesh(), triangulate_, cellSelection(true));
if (debug)
{
......
......@@ -29,6 +29,34 @@ Description
The cuttingPlane algorithm 'cuts' the mesh.
The plane is triangulated by default.
Example of function object partial specification:
\verbatim
surfaces
(
surface1
{
type plane;
planeType pointAndNormal;
pointAndNormalDict
{
...
}
}
);
\endverbatim
Where the sub-entries comprise:
\table
Property | Description | Required | Default
type | plane | yes |
planeType | plane description (pointAndNormal etc) | yes |
triangulate | triangulate faces | no | true
bounds | limit with bounding box | no |
zone | limit to cell zone (name or regex) | no |
zones | limit to cell zones (names, regexs) | no |
coordinateSystem | define plane within given coordinate system | no |
\endtable
Note
Does not actually cut until update() called.
......@@ -66,10 +94,10 @@ class surfMeshSamplePlane
//- The zone or zones in which cutting is to occur
wordRes zoneNames_;
//- Optional bounding box to trim triangles against
//- Optional bounding box to trim against
const boundBox bounds_;
//- Triangulated faces or keep faces as is
//- Triangulate faces or not
const bool triangulate_;
//- Track if the surface needs an update
......@@ -78,18 +106,11 @@ class surfMeshSamplePlane
// Private Member Functions
//- Check and warn if bounding box does not intersect mesh or plane
void checkBoundsIntersection
(
const plane& pln,
const boundBox& meshBb
) const;
//- Define cell selection from zones and bounding box.
// Optionally check and warn if the plane does not intersect
// with the bounds of the mesh (or submesh) or if the bounding box
// does not overlap with the mesh (or submesh)
bitSet cellSelection(const bool warnIntersect=false) const;
bitSet cellSelection(const bool warn=false) const;
//- Sample field on surface
......
......@@ -46,6 +46,9 @@ SourceFiles
namespace Foam
{
// Forward declarations
class polyMesh;
/*---------------------------------------------------------------------------*\
Class cuttingPlane Declaration
\*---------------------------------------------------------------------------*/
......@@ -97,6 +100,35 @@ protected:
);
//- Check and warn if bounding boxes do not intersect,
//- and if the plane does not intersect the bounding boxes
void checkOverlap
(
const word callerName,
const boundBox& meshBounds,
const boundBox& userBounds
) const;
//- Define cell selection from bounding-box and zones.
//
// \param userBounds Optionally user-specified bounding box
// \param zoneNames Optionally user-specified zone names
// \param warn Check and warn if the plane does not
// intersect with the bounds of the mesh (or submesh) or
// if the bounding box does not overlap with the mesh (or submesh)
//
// \return A set of nCells size with the selected cells or an empty
// set if no bounding-box or zones were specified.
bitSet cellSelection
(
const polyMesh& mesh,
const boundBox& userBounds,
const wordRes& zoneNames,
const word callerName,
const bool warn
) const;
public:
// Constructors
......@@ -134,12 +166,18 @@ public:
// Member Functions
//- Return the plane used
//- The plane used.
const plane& planeDesc() const
{
return *this;
}
//- The plane used. Non-const access.
plane& planeDesc()
{
return *this;
}
// Member Operators
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2018 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "cuttingPlane.H"
#include "polyMesh.H"
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void Foam::cuttingPlane::checkOverlap
(
const word callerName,
const boundBox& meshBounds,
const boundBox& userBounds
) const
{
cuttingSurfaceBase::checkOverlap(callerName, meshBounds, userBounds);