Commit a8ef9e97 authored by Mark OLESEN's avatar Mark OLESEN
Browse files

ENH: make cuttingPlane cell walker an algorithm

- takes two general actions:
  1. orient edge in canonical direction (positive gradient) and detect
     any edge intersection.
  2. edge intersection alpha (0-1)

- refactor into a cuttingSurfaceBase intermediate class with the
  actions as templated parameters rather than function pointers. This
  allows the use of lambda functions with captures from the caller.
parent 7cf232ce
......@@ -21,7 +21,7 @@ sampledSet/shortestPath/shortestPathSet.C
surface/cutting/cuttingPlane.C
surface/cutting/cuttingPlaneCuts.C
surface/cutting/cuttingPlaneWalk.C
surface/cutting/cuttingSurfaceBase.C
surface/distanceSurface/distanceSurface.C
surface/isoSurface/isoSurface.C
surface/isoSurface/isoSurfaceCell.C
......
......@@ -25,11 +25,6 @@ License
#include "cuttingPlane.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
int Foam::cuttingPlane::debug(Foam::debug::debugSwitch("cuttingPlane", 0));
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::cuttingPlane::cuttingPlane(const plane& pln)
......@@ -89,6 +84,9 @@ void Foam::cuttingPlane::performCut
bitSet&& cellIdLabels
)
{
const plane& pln = *this;
const pointField& pts = mesh.points();
MeshStorage::clear();
meshCells_.clear();
......@@ -107,57 +105,40 @@ void Foam::cuttingPlane::performCut
// - 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));
}
// Walk cell cuts to create faces
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)
// Action #1:
// - Orient edge so it points in the positive normal direction.
// - Edge/plane intersection when the sign changes
const auto edgeOrientIntersect =
[=](edge& e) -> bool
{
remappedCells[facei] = meshCells_[faceMap[facei]];
}
meshCells_.transfer(remappedCells);
}
if (sides[e.last()] < sides[e.first()])
{
e.flip();
}
return sides[e.first()] != sides[e.last()];
};
// Action #2:
// - The edge intersection alpha
const auto edgeAlphaIntersect =
[=](const edge& e) -> scalar
{
return pln.lineIntersect(e.line(pts));
};
walkCellCuts
(
mesh,
cellCuts,
edgeOrientIntersect,
edgeAlphaIntersect,
triangulate,
nFaceCuts
);
}
......
......@@ -39,18 +39,13 @@ SourceFiles
#define cuttingPlane_H
#include "plane.H"
#include "faceList.H"
#include "bitSet.H"
#include "MeshedSurface.H"
#include "cuttingSurfaceBase.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declarations
class primitiveMesh;
/*---------------------------------------------------------------------------*\
Class cuttingPlane Declaration
\*---------------------------------------------------------------------------*/
......@@ -58,18 +53,8 @@ class primitiveMesh;
class cuttingPlane
:
public plane,
public MeshedSurface<face>
public cuttingSurfaceBase
{
//- Private typedef for convenience
typedef MeshedSurface<face> MeshStorage;
// Private data
//- List of cells cut by the plane
labelList meshCells_;
// Private Member Functions
//- Determine cut cells, possibly restricted to a list of cells
......@@ -88,21 +73,6 @@ class cuttingPlane
bitSet& cellCuts
);
//- Walk the cell cuts to create faces
//
// \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
(
const primitiveMesh& mesh,
const bitSet& cellCuts,
const PackedList<2>& sides,
const bool triangulate,
const label nFaceCuts = 0
);
protected:
......@@ -114,41 +84,21 @@ protected:
// Protected Member Functions
//- Cut mesh with existing plane, restricted to a list of cells
//- Cut mesh, restricted to a list of cells
using cuttingSurfaceBase::performCut;
//- Cut mesh, restricted to a list of cells
// Reclaim memory for cellIdLabels
void performCut
virtual void performCut
(
const primitiveMesh& mesh,
const bool triangulate,
bitSet&& cellIdLabels
);
//- Cut mesh with existing plane, restricted to a list of cells
void performCut
(
const primitiveMesh& mesh,
const bool triangulate,
const bitSet& cellIdLabels = bitSet()
);
//- Cut mesh with existing plane, restricted to a list of cells
void performCut
(
const primitiveMesh& mesh,
const bool triangulate,
const labelUList& cellIdLabels
);
//- Remap action on triangulation or cleanup
virtual void remapFaces(const labelUList& faceMap);
public:
//- Debug information
static int debug;
// Constructors
//- Construct from plane and mesh reference,
......@@ -190,24 +140,6 @@ public:
return *this;
}
//- The mesh cells cut by the plane
const labelList& meshCells() const
{
return meshCells_;
}
//- The mesh cells cut by the plane
labelList& meshCells()
{
return meshCells_;
}
//- Have any cells been cut?
bool cut() const
{
return meshCells_.size();
}
// 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 "cuttingSurfaceBase.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
int Foam::cuttingSurfaceBase::debug
(
Foam::debug::debugSwitch("cuttingSurfaceBase", 0)
);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::cuttingSurfaceBase::cuttingSurfaceBase()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::cuttingSurfaceBase::performCut
(
const primitiveMesh& mesh,
const bool triangulate,
const bitSet& cellIdLabels
)
{
bitSet subsetCells(cellIdLabels);
performCut(mesh, triangulate, std::move(subsetCells));
}
void Foam::cuttingSurfaceBase::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::cuttingSurfaceBase::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::cuttingSurfaceBase::operator=(const cuttingSurfaceBase& rhs)
{
if (this == &rhs)
{
FatalErrorInFunction
<< "Attempted assignment to self"
<< abort(FatalError);
}
static_cast<MeshStorage&>(*this) = rhs;
meshCells_ = rhs.meshCells();
}
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / 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/>.
Class
Foam::cuttingSurfaceBase
Description
Base for creating a MeshedSurface by performing some type of cell
cutting/intersection.
No attempt at resolving degenerate cases.
Since the cut faces can be quite ugly, they will often be triangulated.
SourceFiles
cuttingSurfaceBase.C
\*---------------------------------------------------------------------------*/
#ifndef cuttingSurfaceBase_H
#define cuttingSurfaceBase_H
#include "bitSet.H"
#include "faceList.H"
#include "MeshedSurface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declarations
class primitiveMesh;
/*---------------------------------------------------------------------------*\
Class cuttingSurfaceBase Declaration
\*---------------------------------------------------------------------------*/
class cuttingSurfaceBase
:
public MeshedSurface<face>
{
protected:
//- Typedef for convenience
typedef MeshedSurface<face> MeshStorage;
// Protected Data
//- List of the cells cut
labelList meshCells_;
// Protected Member Functions
//- Walk cell cuts to create faces
//
// \tparam EdgeOrientIntersect
// Parameter (edge&), returns bool.
// Orient edge for a consistent positive gradient.
// Checks for edge intersection (true|false).
//
// \tparam EdgeAlphaIntersect
// Parameter (const edge&), returns scalar.
// Determine alpha [0-1] for an intersecting edge.
// No guarantees when used with non-intersecting edges.
//
// \param cellCuts [in] The cells to walk.
template<class EdgeOrientIntersect, class EdgeAlphaIntersect>
void walkCellCuts
(
const primitiveMesh& mesh,
const bitSet& cellCuts,
const EdgeOrientIntersect& edgeOrientIntersect,
const EdgeAlphaIntersect& edgeAlphaIntersect,
const bool triangulate,
label nFaceCuts = 0
);
//- Cut mesh, restricted to a list of cells
virtual void performCut
(
const primitiveMesh& mesh,
const bool triangulate,
const labelUList& cellIdLabels
);
//- Cut mesh, restricted to a list of cells
virtual void performCut
(
const primitiveMesh& mesh,
const bool triangulate,
const bitSet& cellSelectionMask = bitSet()
);
//- Cut mesh, restricted to a list of cells
// Reclaim memory for cellSelectionMask
virtual void performCut
(
const primitiveMesh& mesh,
const bool triangulate,
bitSet&& cellSelectionMask
) = 0;
//- Remap action on triangulation or cleanup
virtual void remapFaces(const labelUList& faceMap);
public:
//- Debug information
static int debug;
// Constructors
//- Construct null
cuttingSurfaceBase();
//- Destructors
virtual ~cuttingSurfaceBase() = default;
// Member Functions
//- The mesh cells cut
const labelList& meshCells() const
{
return meshCells_;
}
//- The mesh cells cut
labelList& meshCells()
{
return meshCells_;
}
//- Have any cells been cut?
bool cut() const
{
return meshCells_.size();
}
// Member Operators
//- Copy assignment
void operator=(const cuttingSurfaceBase& rhs);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "cuttingSurfaceBaseTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
......@@ -23,40 +23,19 @@ License
\*---------------------------------------------------------------------------*/
#include "cuttingPlane.H"
#include "volFields.H"
#include "edgeHashes.H"
#include "HashOps.H"
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
namespace Foam
{
// Check edge/plane intersection based on crossings ... trivial check.
// Orients the edge (first,last) points in the positive normal direction
inline bool intersectEdgeOrient(const PackedList<2>& sides, edge& e)
{
if (sides[e.first()] == sides[e.last()])
{
return false;
}
else if (sides[e.last()] < sides[e.first()])
{
e.flip();
}
return true;
}
} // End namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::cuttingPlane::walkCellCuts
template<class EdgeOrientIntersect, class EdgeAlphaIntersect>
void Foam::cuttingSurfaceBase::walkCellCuts
(
const primitiveMesh& mesh,
const bitSet& cellCuts,
const PackedList<2>& sides,
const EdgeOrientIntersect& edgeOrientIntersect,
const EdgeAlphaIntersect& edgeAlphaIntersect,
const bool triangulate,
label nFaceCuts
)
......@@ -66,7 +45,6 @@ void Foam::cuttingPlane::walkCellCuts
const cellList& cells = mesh.cells();
const pointField& points = mesh.points();
// Dynamic lists to handle triangulation and/or missed cuts etc
const label nCellCuts = cellCuts.count();
......@@ -156,8 +134,8 @@ void Foam::cuttingPlane::walkCellCuts
{
edge e(f.faceEdge(fp));
// Action #1: orient edge and detect intersection
if (!intersectEdgeOrient(sides, e))
// Action #1: Orient edge (+ve gradient) and detect intersect
if (!edgeOrientIntersect(e))
{
continue;
}
......@@ -191,8 +169,7 @@ void Foam::cuttingPlane::walkCellCuts
const point& p1 = points[e.last()];
// Action #2: edge cut alpha
const scalar alpha =
this->lineIntersect(linePointRef(p0, p1));
const scalar alpha = edgeAlphaIntersect(e);
if (alpha < SMALL)