STYLE: split up cuttingPlane source files

This commit is contained in:
Mark Olesen
2018-09-14 09:19:32 +02:00
parent 53b0bb0782
commit 7cf232ceec
5 changed files with 411 additions and 349 deletions

View File

@ -19,7 +19,9 @@ sampledSet/uniform/uniformSet.C
sampledSet/array/arraySet.C sampledSet/array/arraySet.C
sampledSet/shortestPath/shortestPathSet.C sampledSet/shortestPath/shortestPathSet.C
surface/cuttingPlane/cuttingPlane.C surface/cutting/cuttingPlane.C
surface/cutting/cuttingPlaneCuts.C
surface/cutting/cuttingPlaneWalk.C
surface/distanceSurface/distanceSurface.C surface/distanceSurface/distanceSurface.C
surface/isoSurface/isoSurface.C surface/isoSurface/isoSurface.C
surface/isoSurface/isoSurfaceCell.C surface/isoSurface/isoSurfaceCell.C

View File

@ -0,0 +1,181 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
int Foam::cuttingPlane::debug(Foam::debug::debugSwitch("cuttingPlane", 0));
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::cuttingPlane::cuttingPlane(const plane& pln)
:
plane(pln)
{}
Foam::cuttingPlane::cuttingPlane
(
const plane& pln,
const primitiveMesh& mesh,
const bool triangulate,
const bitSet& cellIdLabels
)
:
plane(pln)
{
performCut(mesh, triangulate, cellIdLabels);
}
Foam::cuttingPlane::cuttingPlane
(
const plane& pln,
const primitiveMesh& mesh,
const bool triangulate,
bitSet&& cellIdLabels
)
:
plane(pln)
{
performCut(mesh, triangulate, cellIdLabels);
}
Foam::cuttingPlane::cuttingPlane
(
const plane& pln,
const primitiveMesh& mesh,
const bool triangulate,
const labelUList& cellIdLabels
)
:
plane(pln)
{
performCut(mesh, triangulate, cellIdLabels);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::cuttingPlane::performCut
(
const primitiveMesh& mesh,
const bool triangulate,
bitSet&& cellIdLabels
)
{
MeshStorage::clear();
meshCells_.clear();
// Pre-populate with restriction
bitSet cellCuts(std::move(cellIdLabels));
if (cellCuts.size())
{
cellCuts.resize(mesh.nCells());
}
// For each mesh point, the encoded side (0,1,2) of the plane
PackedList<2> sides;
// Determine cells that are (likely) cut
// - some ambiguity when plane is exactly between cells
const label nFaceCuts = calcCellCuts(mesh, sides, cellCuts);
// Find closed loop from cell cuts
walkCellCuts(mesh, cellCuts, sides, triangulate, nFaceCuts);
}
void Foam::cuttingPlane::performCut
(
const primitiveMesh& mesh,
const bool triangulate,
const bitSet& cellIdLabels
)
{
bitSet subsetCells(cellIdLabels);
performCut(mesh, triangulate, std::move(subsetCells));
}
void Foam::cuttingPlane::performCut
(
const primitiveMesh& mesh,
const bool triangulate,
const labelUList& cellIdLabels
)
{
bitSet subsetCells;
if (notNull(cellIdLabels))
{
// Pre-populate with restriction
subsetCells.resize(mesh.nCells());
subsetCells.set(cellIdLabels);
}
performCut(mesh, triangulate, std::move(subsetCells));
}
void Foam::cuttingPlane::remapFaces(const labelUList& faceMap)
{
if (notNull(faceMap) && !faceMap.empty())
{
MeshStorage::remapFaces(faceMap);
List<label> remappedCells(faceMap.size());
forAll(faceMap, facei)
{
remappedCells[facei] = meshCells_[faceMap[facei]];
}
meshCells_.transfer(remappedCells);
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void Foam::cuttingPlane::operator=(const cuttingPlane& rhs)
{
if (this == &rhs)
{
FatalErrorInFunction
<< "Attempted assignment to self"
<< abort(FatalError);
}
static_cast<MeshStorage&>(*this) = rhs;
static_cast<plane&>(*this) = rhs;
meshCells_ = rhs.meshCells();
}
// ************************************************************************* //

View File

@ -30,10 +30,6 @@ Description
No attempt at resolving degenerate cases. No attempt at resolving degenerate cases.
Since the cut faces can be quite ugly, they will often be triangulated. Since the cut faces can be quite ugly, they will often be triangulated.
Note
When the cutting plane coincides with a mesh face, the cell edge on the
positive side of the plane is taken.
SourceFiles SourceFiles
cuttingPlane.C cuttingPlane.C
@ -43,7 +39,6 @@ SourceFiles
#define cuttingPlane_H #define cuttingPlane_H
#include "plane.H" #include "plane.H"
#include "pointField.H"
#include "faceList.H" #include "faceList.H"
#include "bitSet.H" #include "bitSet.H"
#include "MeshedSurface.H" #include "MeshedSurface.H"
@ -79,26 +74,31 @@ class cuttingPlane
//- Determine cut cells, possibly restricted to a list of cells //- Determine cut cells, possibly restricted to a list of cells
// //
// \param sides [out] For each mesh point, the encoded side of the
// plane (0=BACK, 1=ONPLANE, 2=FRONT).
// \param cellCuts [in,out] On input an empty set (ie, no restriction) // \param cellCuts [in,out] On input an empty set (ie, no restriction)
// or subsetted cells. On output, the cells cut according to the // or subsetted cells. On output, the cells cut according to the
// planeSides detection. // plane-sides detection.
// //
// \return number of faces cut // \return number of faces cut
label calcCellCuts label calcCellCuts
( (
const primitiveMesh& mesh, const primitiveMesh& mesh,
const PackedList<2>& planeSides, PackedList<2>& sides,
bitSet& cellCuts bitSet& cellCuts
); );
//- Walk the cell cuts to create faces //- Walk the cell cuts to create faces
// //
// \param planeSides [in] Used to determine edge cuts // \param cellCuts [in] The cells to walk.
// \param sides [in] For each mesh point, the encoded side of the
// plane (0=BACK, 1=ONPLANE, 2=FRONT), which is used to determine
// the edge cuts
void walkCellCuts void walkCellCuts
( (
const primitiveMesh& mesh, const primitiveMesh& mesh,
const bitSet& cellCuts, const bitSet& cellCuts,
const PackedList<2>& planeSides, const PackedList<2>& sides,
const bool triangulate, const bool triangulate,
const label nFaceCuts = 0 const label nFaceCuts = 0
); );
@ -115,20 +115,20 @@ protected:
// Protected Member Functions // Protected Member Functions
//- Cut mesh with existing plane, restricted to a list of cells //- Cut mesh with existing plane, restricted to a list of cells
// Reclaim memory for cellIdLabels
void performCut void performCut
( (
const primitiveMesh& mesh, const primitiveMesh& mesh,
const bool triangulate, const bool triangulate,
const bitSet& cellSelectionMask = bitSet() bitSet&& cellIdLabels
); );
//- Cut mesh with existing plane, restricted to a list of cells //- Cut mesh with existing plane, restricted to a list of cells
// Reclaim memory for cellSelectionMask
void performCut void performCut
( (
const primitiveMesh& mesh, const primitiveMesh& mesh,
const bool triangulate, const bool triangulate,
bitSet&& cellSelectionMask const bitSet& cellIdLabels = bitSet()
); );
//- Cut mesh with existing plane, restricted to a list of cells //- Cut mesh with existing plane, restricted to a list of cells

View File

@ -0,0 +1,196 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "volFields.H"
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
namespace Foam
{
// Check for face/plane intersection based on crossings
// Took (-1,0,+1) from plane::sign and packed as (0,1,2).
// Now use for left shift to obtain (1,2,4).
//
// Test accumulated value for an intersection with the plane.
inline bool intersectsFace
(
const PackedList<2>& sides,
const labelUList& indices
)
{
unsigned accum = 0u;
for (const label pointi : indices)
{
accum |= (1u << sides[pointi]);
}
// Accumulated value 3,5,6,7 indicates an intersection
return (accum == 3 || accum >= 5);
}
} // End namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::label Foam::cuttingPlane::calcCellCuts
(
const primitiveMesh& mesh,
PackedList<2>& sides,
bitSet& cellCuts
)
{
const faceList& faces = mesh.faces();
const pointField& pts = mesh.points();
const label nCells = mesh.nCells();
const label nFaces = mesh.nFaces();
const label nInternFaces = mesh.nInternalFaces();
// Classify sides of plane (0=BACK, 1=ONPLANE, 2=FRONT) for each point
{
const plane& pln = *this;
const label len = pts.size();
sides.resize(len);
// From signed (-1,0,+1) to (0,1,2) for PackedList
for (label i=0; i < len; ++i)
{
sides.set(i, unsigned(1 + pln.sign(pts[i], SMALL)));
}
}
// Detect cells cuts from the face cuts
label nFaceCuts = 0;
// 1st face cut of cell
bitSet hasCut1(nCells);
// 2nd face cut of cell
bitSet hasCut2(nCells);
for (label facei = 0; facei < nInternFaces; ++facei)
{
if (intersectsFace(sides, faces[facei]))
{
const label own = mesh.faceOwner()[facei];
const label nei = mesh.faceNeighbour()[facei];
++nFaceCuts;
if (!hasCut1.set(own))
{
hasCut2.set(own);
}
if (!hasCut1.set(nei))
{
hasCut2.set(nei);
}
}
}
for (label facei = nInternFaces; facei < nFaces; ++facei)
{
if (intersectsFace(sides, faces[facei]))
{
const label own = mesh.faceOwner()[facei];
++nFaceCuts;
if (!hasCut1.set(own))
{
hasCut2.set(own);
}
}
}
hasCut1.clearStorage(); // Not needed now
if (cellCuts.size())
{
cellCuts.resize(nCells); // safety
cellCuts &= hasCut2; // restrict to cell subset
if (debug)
{
Pout<<"detected " << cellCuts.count() << "/" << nCells
<< " cells cut, subsetted from "
<< hasCut2.count() << "/" << nCells << " cells." << endl;
}
}
else
{
cellCuts = std::move(hasCut2);
if (debug)
{
Pout<<"detected " << cellCuts.count() << "/" << nCells
<< " cells cut." << endl;
}
}
if (debug && isA<fvMesh>(mesh))
{
const fvMesh& fvm = dynamicCast<const fvMesh&>(mesh);
volScalarField cCuts
(
IOobject
(
"cuttingPlane.cellCuts",
fvm.time().timeName(),
fvm.time(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
fvm,
dimensionedScalar(dimless, Zero)
);
auto& cCutsFld = cCuts.primitiveFieldRef();
for (const label celli : cellCuts)
{
cCutsFld[celli] = 1;
}
Pout<< "Writing cut types:" << cCuts.objectPath() << endl;
cCuts.write();
}
return nFaceCuts;
}
// ************************************************************************* //

View File

@ -2,8 +2,8 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2018 OpenCFD Ltd.
\\/ M anipulation | Copyright (C) 2018 OpenCFD Ltd. \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
This file is part of OpenFOAM. This file is part of OpenFOAM.
@ -24,17 +24,9 @@ License
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "cuttingPlane.H" #include "cuttingPlane.H"
#include "fvMesh.H"
#include "volFields.H"
#include "meshTools.H"
#include "edgeHashes.H" #include "edgeHashes.H"
#include "HashOps.H" #include "HashOps.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
int Foam::cuttingPlane::debug(Foam::debug::debugSwitch("cuttingPlane", 0));
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
namespace Foam namespace Foam
@ -55,172 +47,11 @@ namespace Foam
return true; return true;
} }
// Classify sides of plane (0=BACK, 1=ONPLANE, 2=FRONT) for each point
inline PackedList<2> classifySides(const plane& pln, const pointField& pts)
{
const label len = pts.size();
PackedList<2> output(len);
// From signed (-1,0,+1) to (0,1,2) for PackedList
for (label i=0; i < len; ++i)
{
output.set(i, unsigned(1 + pln.sign(pts[i], SMALL)));
}
return output;
}
// Check for face/plane intersection based on crossings
// Took (-1,0,+1) from plane::sign and packed as (0,1,2).
// Now use for left shift to obtain (1,2,4).
//
// Test accumulated value for an intersection with the plane.
inline bool intersectsFace
(
const PackedList<2>& sides,
const labelUList& indices
)
{
unsigned accum = 0u;
for (const label pointi : indices)
{
accum |= (1u << sides[pointi]);
}
// Accumulated value 3,5,6,7 indicates an intersection
return (accum == 3 || accum >= 5);
}
} // End namespace Foam } // End namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::label Foam::cuttingPlane::calcCellCuts
(
const primitiveMesh& mesh,
const PackedList<2>& sides,
bitSet& cellCuts
)
{
const faceList& faces = mesh.faces();
const label nCells = mesh.nCells();
const label nFaces = mesh.nFaces();
const label nInternFaces = mesh.nInternalFaces();
// Detect cells cuts from the face cuts
label nFaceCuts = 0;
// 1st face cut of cell
bitSet hasCut1(nCells);
// 2nd face cut of cell
bitSet hasCut2(nCells);
for (label facei = 0; facei < nInternFaces; ++facei)
{
if (intersectsFace(sides, faces[facei]))
{
const label own = mesh.faceOwner()[facei];
const label nei = mesh.faceNeighbour()[facei];
++nFaceCuts;
if (!hasCut1.set(own))
{
hasCut2.set(own);
}
if (!hasCut1.set(nei))
{
hasCut2.set(nei);
}
}
}
for (label facei = nInternFaces; facei < nFaces; ++facei)
{
if (intersectsFace(sides, faces[facei]))
{
const label own = mesh.faceOwner()[facei];
++nFaceCuts;
if (!hasCut1.set(own))
{
hasCut2.set(own);
}
}
}
hasCut1.clearStorage(); // Not needed now
if (cellCuts.size())
{
cellCuts.resize(nCells); // safety
cellCuts &= hasCut2; // restrict to cell subset
if (debug)
{
Pout<<"detected " << cellCuts.count() << "/" << nCells
<< " cells cut, subsetted from "
<< hasCut2.count() << "/" << nCells << " cells." << endl;
}
}
else
{
cellCuts = std::move(hasCut2);
if (debug)
{
Pout<<"detected " << cellCuts.count() << "/" << nCells
<< " cells cut." << endl;
}
}
if (debug && isA<fvMesh>(mesh))
{
const fvMesh& fvm = dynamicCast<const fvMesh&>(mesh);
volScalarField debugField
(
IOobject
(
"cuttingPlane.cellCuts",
fvm.time().timeName(),
fvm.time(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
fvm,
dimensionedScalar(dimless, Zero)
);
auto& debugFld = debugField.primitiveFieldRef();
for (const label celli : cellCuts)
{
debugFld[celli] = 1;
}
Pout<< "Writing cut types:"
<< debugField.objectPath() << endl;
debugField.write();
}
return nFaceCuts;
}
void Foam::cuttingPlane::walkCellCuts void Foam::cuttingPlane::walkCellCuts
( (
const primitiveMesh& mesh, const primitiveMesh& mesh,
@ -325,7 +156,7 @@ void Foam::cuttingPlane::walkCellCuts
{ {
edge e(f.faceEdge(fp)); edge e(f.faceEdge(fp));
// Action #1: detect edge intersection and orient edge // Action #1: orient edge and detect intersection
if (!intersectEdgeOrient(sides, e)) if (!intersectEdgeOrient(sides, e))
{ {
continue; continue;
@ -356,8 +187,8 @@ void Foam::cuttingPlane::walkCellCuts
// Expected id for the cut point // Expected id for the cut point
cutPointId = dynCutPoints.size(); cutPointId = dynCutPoints.size();
const point& p0 = points[e[0]]; const point& p0 = points[e.first()];
const point& p1 = points[e[1]]; const point& p1 = points[e.last()];
// Action #2: edge cut alpha // Action #2: edge cut alpha
const scalar alpha = const scalar alpha =
@ -365,41 +196,41 @@ void Foam::cuttingPlane::walkCellCuts
if (alpha < SMALL) if (alpha < SMALL)
{ {
// -> equal(alpha, 0) with more tolerance pointCutType |= 0x1; // Cut at 0 (first)
if (endPoints.insert(e[0], cutPointId)) const label endp = e.first();
if (endPoints.insert(endp, cutPointId))
{ {
localEndPoints.insert(e[0]); localEndPoints.insert(endp);
dynCutPoints.append(p0); dynCutPoints.append(p0);
} }
else else
{ {
cutPointId = endPoints[e[0]]; cutPointId = endPoints[endp];
} }
pointCutType |= 0x1; // Cut at 0
} }
else if (alpha >= (1.0 - SMALL)) else if (alpha >= (1.0 - SMALL))
{ {
// -> equal(alpha, 1) with more tolerance pointCutType |= 0x2; // Cut at 1 (last)
if (endPoints.insert(e[1], cutPointId)) const label endp = e.last();
if (endPoints.insert(endp, cutPointId))
{ {
localEndPoints.insert(e[1]); localEndPoints.insert(endp);
dynCutPoints.append(p1); dynCutPoints.append(p1);
} }
else else
{ {
cutPointId = endPoints[e[1]]; cutPointId = endPoints[endp];
} }
pointCutType |= 0x2; // Cut at 1
} }
else else
{ {
dynCutPoints.append((1-alpha)*p0 + alpha*p1);
pointCutType |= 0x4; // Cut between pointCutType |= 0x4; // Cut between
dynCutPoints.append((1-alpha)*p0 + alpha*p1);
} }
// Introduce new edge cut point // Introduce new edge cut point
@ -565,152 +396,4 @@ void Foam::cuttingPlane::walkCellCuts
} }
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::cuttingPlane::cuttingPlane(const plane& pln)
:
plane(pln)
{}
Foam::cuttingPlane::cuttingPlane
(
const plane& pln,
const primitiveMesh& mesh,
const bool triangulate,
const bitSet& cellIdLabels
)
:
plane(pln)
{
performCut(mesh, triangulate, cellIdLabels);
}
Foam::cuttingPlane::cuttingPlane
(
const plane& pln,
const primitiveMesh& mesh,
const bool triangulate,
bitSet&& cellIdLabels
)
:
plane(pln)
{
performCut(mesh, triangulate, cellIdLabels);
}
Foam::cuttingPlane::cuttingPlane
(
const plane& pln,
const primitiveMesh& mesh,
const bool triangulate,
const labelUList& cellIdLabels
)
:
plane(pln)
{
performCut(mesh, triangulate, cellIdLabels);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::cuttingPlane::performCut
(
const primitiveMesh& mesh,
const bool triangulate,
bitSet&& cellIdLabels
)
{
MeshStorage::clear();
meshCells_.clear();
// Pre-populate with restriction
bitSet cellCuts(std::move(cellIdLabels));
if (cellCuts.size())
{
cellCuts.resize(mesh.nCells());
}
// Classification of points with respect to the plane
PackedList<2> sides = classifySides(*this, mesh.points());
// Determine cells that are (likely) cut
// - some ambiguity when plane is exactly between cells
const label nFaceCuts = calcCellCuts(mesh, sides, cellCuts);
// Find closed loop from cell cuts
walkCellCuts(mesh, cellCuts, sides, triangulate, nFaceCuts);
}
void Foam::cuttingPlane::performCut
(
const primitiveMesh& mesh,
const bool triangulate,
const bitSet& cellIdLabels
)
{
bitSet subsetCells(cellIdLabels);
performCut(mesh, triangulate, std::move(subsetCells));
}
void Foam::cuttingPlane::performCut
(
const primitiveMesh& mesh,
const bool triangulate,
const labelUList& cellIdLabels
)
{
bitSet subsetCells;
if (notNull(cellIdLabels))
{
// Pre-populate with restriction
subsetCells.resize(mesh.nCells());
subsetCells.set(cellIdLabels);
}
performCut(mesh, triangulate, std::move(subsetCells));
}
void Foam::cuttingPlane::remapFaces(const labelUList& faceMap)
{
if (notNull(faceMap) && !faceMap.empty())
{
MeshStorage::remapFaces(faceMap);
List<label> newCutCells(faceMap.size());
forAll(faceMap, facei)
{
newCutCells[facei] = meshCells_[faceMap[facei]];
}
meshCells_.transfer(newCutCells);
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void Foam::cuttingPlane::operator=(const cuttingPlane& rhs)
{
if (this == &rhs)
{
FatalErrorInFunction
<< "Attempted assignment to self"
<< abort(FatalError);
}
static_cast<MeshStorage&>(*this) = rhs;
static_cast<plane&>(*this) = rhs;
meshCells_ = rhs.meshCells();
}
// ************************************************************************* // // ************************************************************************* //