Commit de2eed3e authored by Mark Olesen's avatar Mark Olesen
Browse files

ENH: additional methods and improvements to plane

- signedDistance() method is like distance() but retains
  the positive/negative sign for the side of the plane.

- the sign() method returns the sign as -1,0,+1 integer for
  classification purposes where it is important to distinguish between
  a zero value and a positive value (eg, for cutting). Optional
  tolerance can be supplied to round for zero.

- refactor and inlined simple and frequently used methods.

- add boundBox faceCentre() method, which can be useful for creating
  clipping planes from a bounding box.
  Relocated treeBoundBox faceNormals to boundBox since they apply
  equally there - the meaning of the faces (x-min, x-max, etc)
  is the same, even if the point addressing for the faces differs.
parent 84e2df49
......@@ -99,11 +99,11 @@ bool Foam::conformalVoronoiMesh::meshableRegion
{
case extendedFeatureEdgeMesh::INSIDE:
{
return (side == plane::FLIP) ? true : false;
return (side == plane::BACK);
}
case extendedFeatureEdgeMesh::OUTSIDE:
{
return (side == plane::NORMAL) ? true : false;
return (side == plane::FRONT);
}
case extendedFeatureEdgeMesh::BOTH:
{
......@@ -132,12 +132,12 @@ bool Foam::conformalVoronoiMesh::regionIsInside
{
plane::side sideA
(
((masterPtVec & normalA) <= 0) ? plane::FLIP : plane::NORMAL
((masterPtVec & normalA) <= 0) ? plane::BACK : plane::FRONT
);
plane::side sideB
(
((masterPtVec & normalB) <= 0) ? plane::FLIP : plane::NORMAL
((masterPtVec & normalB) <= 0) ? plane::BACK : plane::FRONT
);
const bool meshableRegionA = meshableRegion(sideA, volTypeA);
......
......@@ -18,8 +18,8 @@ planeType pointAndNormal;
pointAndNormalDict
{
basePoint (0 0 0);
normalVector (0 1 0);
point (0 0 0);
normal (0 1 0);
}
planeTolerance 1e-3;
......
......@@ -493,7 +493,7 @@ int main(int argc, char *argv[])
Info<< "Only include feature edges that intersect the plane"
<< " with normal " << cutPlane.normal()
<< " and base point " << cutPlane.refPoint() << endl;
<< " and origin " << cutPlane.origin() << endl;
features().subsetPlane(edgeStat, cutPlane);
}
......
......@@ -45,12 +45,22 @@ const Foam::boundBox Foam::boundBox::invertedBox
const Foam::faceList Foam::boundBox::faces
({
// Point and face order as per hex cellmodel
face{0, 4, 7, 3}, // x-min
face{1, 2, 6, 5}, // x-max
face{0, 1, 5, 4}, // y-min
face{3, 7, 6, 2}, // y-max
face{0, 3, 2, 1}, // z-min
face{4, 5, 6, 7} // z-max
face({0, 4, 7, 3}), // 0: x-min, left
face({1, 2, 6, 5}), // 1: x-max, right
face({0, 1, 5, 4}), // 2: y-min, bottom
face({3, 7, 6, 2}), // 3: y-max, top
face({0, 3, 2, 1}), // 4: z-min, back
face({4, 5, 6, 7}) // 5: z-max, front
});
const Foam::FixedList<Foam::vector, 6> Foam::boundBox::faceNormals
({
vector(-1, 0, 0), // 0: x-min, left
vector( 1, 0, 0), // 1: x-max, right
vector( 0, -1, 0), // 2: y-min, bottom
vector( 0, 1, 0), // 3: y-max, top
vector( 0, 0, -1), // 4: z-min, back
vector( 0, 0, 1) // 5: z-max, front
});
......@@ -107,8 +117,8 @@ Foam::boundBox::boundBox
Foam::tmp<Foam::pointField> Foam::boundBox::points() const
{
tmp<pointField> tpoints(new pointField(8));
pointField& pt = tpoints.ref();
auto tpt = tmp<pointField>::New(8);
auto& pt = tpt.ref();
pt[0] = min_; // min-x, min-y, min-z
pt[1] = point(max_.x(), min_.y(), min_.z()); // max-x, min-y, min-z
......@@ -119,7 +129,46 @@ Foam::tmp<Foam::pointField> Foam::boundBox::points() const
pt[6] = max_; // max-x, max-y, max-z
pt[7] = point(min_.x(), max_.y(), max_.z()); // min-x, max-y, max-z
return tpoints;
return tpt;
}
Foam::tmp<Foam::pointField> Foam::boundBox::faceCentres() const
{
auto tpts = tmp<pointField>::New(6);
auto& pts = tpts.ref();
forAll(pts, facei)
{
pts[facei] = faceCentre(facei);
}
return tpts;
}
Foam::point Foam::boundBox::faceCentre(const direction facei) const
{
point pt = boundBox::midpoint();
if (facei > 5)
{
FatalErrorInFunction
<< "face should be [0..5]"
<< abort(FatalError);
}
switch (facei)
{
case 0: pt.x() = min().x(); break; // 0: x-min, left
case 1: pt.x() = max().x(); break; // 1: x-max, right
case 2: pt.y() = min().y(); break; // 2: y-min, bottom
case 3: pt.y() = max().y(); break; // 3: y-max, top
case 4: pt.z() = min().z(); break; // 4: z-min, back
case 5: pt.z() = max().z(); break; // 5: z-max, front
}
return pt;
}
......@@ -160,11 +209,11 @@ bool Foam::boundBox::intersects(const plane& pln) const
bool below = false;
tmp<pointField> tpts(points());
const pointField& pts = tpts();
const auto& pts = tpts();
for (const point& p : pts)
{
if (pln.sideOfPlane(p) == plane::NORMAL)
if (pln.sideOfPlane(p) == plane::FRONT)
{
above = true;
}
......@@ -208,9 +257,9 @@ bool Foam::boundBox::contains
return true;
}
forAll(indices, i)
for (const label pointi : indices)
{
if (!contains(points[indices[i]]))
if (!contains(points[pointi]))
{
return false;
}
......@@ -250,9 +299,9 @@ bool Foam::boundBox::containsAny
return true;
}
forAll(indices, i)
for (const label pointi : indices)
{
if (contains(points[indices[i]]))
if (contains(points[pointi]))
{
return true;
}
......
......@@ -78,6 +78,9 @@ public:
//- Faces to point addressing, as per a 'hex' cell
static const faceList faces;
//- The unit normal per face
static const FixedList<vector, 6> faceNormals;
// Constructors
......@@ -177,6 +180,12 @@ public:
//- Corner points in an order corresponding to a 'hex' cell
tmp<pointField> points() const;
//- Face midpoints
tmp<pointField> faceCentres() const;
//- Face centre of given face index
point faceCentre(const direction facei) const;
// Manipulate
......
......@@ -203,9 +203,9 @@ inline void Foam::boundBox::add
{
if (!points.empty())
{
forAll(indices, i)
for (const label pointi : indices)
{
add(points[indices[i]]);
add(points[pointi]);
}
}
}
......
......@@ -74,9 +74,9 @@ void Foam::boundBox::add
// points may be empty, but a FixedList is never empty
if (!points.empty())
{
for (unsigned i=0; i < Size; ++i)
for (const label pointi : indices)
{
add(points[indices[i]]);
add(points[pointi]);
}
}
}
......@@ -95,9 +95,9 @@ bool Foam::boundBox::contains
return false;
}
forAll(indices, i)
for (const label pointi : indices)
{
if (!contains(points[indices[i]]))
if (!contains(points[pointi]))
{
return false;
}
......@@ -120,9 +120,9 @@ bool Foam::boundBox::containsAny
return false;
}
forAll(indices, i)
for (const label pointi : indices)
{
if (contains(points[indices[i]]))
if (contains(points[pointi]))
{
return true;
}
......
......@@ -57,8 +57,7 @@ SourceFiles
namespace Foam
{
// Forward declaration of friend functions and operators
// Forward declarations
class face;
class triFace;
......
......@@ -21,11 +21,11 @@ License
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Description:
Functions which cut triangles and tetrahedra. Generic operations are
applied to each half of a cut.
Description
Functions for cutting triangles and tetrahedra.
Generic operations are applied to each half of a cut.
SourceFiles:
SourceFiles
cutI.H
cutTemplates.C
......@@ -478,9 +478,9 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//- Cut a triangle along the zero plane defined by the given levels. Templated
// on aboveOp and belowOp; the operations applied to the regions on either side
// of the cut.
//- Cut a triangle along the zero plane defined by the given levels.
// Templated on aboveOp and belowOp; the operations applied to the regions
// on either side of the cut.
template<class AboveOp, class BelowOp>
typename cut::opAddResult<AboveOp, BelowOp>::type triCut
(
......@@ -495,7 +495,7 @@ template<class AboveOp, class BelowOp>
typename cut::opAddResult<AboveOp, BelowOp>::type triCut
(
const FixedList<point, 3>& tri,
const plane& s,
const plane& pln,
const AboveOp& aboveOp,
const BelowOp& belowOp
);
......@@ -515,7 +515,7 @@ template<class AboveOp, class BelowOp>
typename cut::opAddResult<AboveOp, BelowOp>::type tetCut
(
const FixedList<point, 4>& tet,
const plane& s,
const plane& pln,
const AboveOp& aboveOp,
const BelowOp& belowOp
);
......
......@@ -51,7 +51,7 @@ typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::triCut
// opposite the first vertex. This may change the sign of the tri.
FixedList<label, 3> indices({0, 1, 2});
label i;
for (i = 0; i < 3; ++ i)
for (i = 0; i < 3; ++i)
{
if (level[(i + 1)%3]*level[(i + 2)%3] >= 0)
{
......@@ -80,7 +80,7 @@ typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::triCut
// Slice off one corner to form a tri and a quad
Pair<scalar> f;
for (label i = 0; i < 2; ++ i)
for (label i = 0; i < 2; ++i)
{
f[i] = l[0]/(l[0] - l[i+1]);
}
......@@ -99,7 +99,7 @@ template<class AboveOp, class BelowOp>
typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::triCut
(
const FixedList<point, 3>& tri,
const plane& p,
const plane& pln,
const AboveOp& aboveOp,
const BelowOp& belowOp
)
......@@ -108,8 +108,7 @@ typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::triCut
FixedList<scalar, 3> level;
for (label i = 0; i < 3; ++i)
{
// uses unit-normal
level[i] = (tri[i] - p.refPoint()) & p.normal();
level[i] = pln.signedDistance(tri[i]);
}
// Run the level set function
......@@ -129,7 +128,7 @@ typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::tetCut
// Get the min and max over all four vertices and quick return if there is
// no change of sign
scalar levelMin = VGREAT, levelMax = - VGREAT;
for (label i = 0; i < 4; ++ i)
for (label i = 0; i < 4; ++i)
{
levelMin = min(levelMin, level[i]);
levelMax = max(levelMax, level[i]);
......@@ -176,7 +175,7 @@ typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::tetCut
if (n > 2)
{
n = 4 - n;
for (label i = 0; i < 2; ++ i)
for (label i = 0; i < 2; ++i)
{
Swap(indices[i], indices[3-i]);
}
......@@ -199,7 +198,7 @@ typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::tetCut
{
// Slice off one corner to form a tet and a prism
FixedList<scalar, 3> f;
for (label i = 0; i < 3; ++ i)
for (label i = 0; i < 3; ++i)
{
f[i] = l[0]/(l[0] - l[i+1]);
}
......@@ -216,9 +215,9 @@ typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::tetCut
{
// Slice off two corners to form two prisms
FixedList<scalar, 4> f;
for (label i = 0; i < 2; ++ i)
for (label i = 0; i < 2; ++i)
{
for (label j = 0; j < 2; ++ j)
for (label j = 0; j < 2; ++j)
{
f[2*i+j] = l[i]/(l[i] - l[j+2]);
}
......@@ -245,7 +244,7 @@ template<class AboveOp, class BelowOp>
typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::tetCut
(
const FixedList<point, 4>& tet,
const plane& p,
const plane& pln,
const AboveOp& aboveOp,
const BelowOp& belowOp
)
......@@ -254,8 +253,7 @@ typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::tetCut
FixedList<scalar, 4> level;
for (label i = 0; i < 4; ++i)
{
// uses unit-normal
level[i] = (tet[i] - p.refPoint()) & p.normal();
level[i] = pln.signedDistance(tet[i]);
}
// Run the level set function
......
......@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016-2017 OpenCFD Ltd.
\\/ M anipulation | Copyright (C) 2016-2018 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -28,51 +28,72 @@ License
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::plane::calcPntAndVec(const scalarList& C)
void Foam::plane::makeUnitNormal
(
const char * const caller,
const bool notTest
)
{
if (mag(C[0]) > VSMALL)
const scalar magNormal(Foam::mag(normal_));
if (magNormal < VSMALL)
{
point_ = vector((-C[3]/C[0]), 0, 0);
FatalErrorInFunction
<< "Plane normal has zero length.\nCalled from " << caller
<< abort(FatalError);
}
else if (mag(C[1]) > VSMALL)
if (notTest)
{
point_ = vector(0, (-C[3]/C[1]), 0);
normal_ /= magNormal;
}
else if (mag(C[2]) > VSMALL)
}
void Foam::plane::calcFromCoeffs
(
const scalar a,
const scalar b,
const scalar c,
const scalar d,
const char * const caller
)
{
if (mag(a) > VSMALL)
{
point_ = vector(0, 0, (-C[3]/C[2]));
origin_ = vector((-d/a), 0, 0);
}
else
else if (mag(b) > VSMALL)
{
FatalErrorInFunction
<< "At least one plane coefficient must have a value"
<< abort(FatalError);
origin_ = vector(0, (-d/b), 0);
}
normal_ = vector(C[0], C[1], C[2]);
scalar magUnitVector(mag(normal_));
if (magUnitVector < VSMALL)
else if (mag(c) > VSMALL)
{
origin_ = vector(0, 0, (-d/c));
}
else
{
FatalErrorInFunction
<< "Plane normal defined with zero length"
<< "At least one plane coefficient must have a value"
<< abort(FatalError);
}
normal_ /= magUnitVector;
normal_ = vector(a, b, c);
makeUnitNormal(caller);
}
void Foam::plane::calcPntAndVec
void Foam::plane::calcFromEmbeddedPoints
(
const point& point1,
const point& point2,
const point& point3
const point& point3,
const char * const caller
)
{
point_ = (point1 + point2 + point3)/3;
vector line12 = point1 - point2;
vector line23 = point2 - point3;
origin_ = (point1 + point2 + point3)/3;
const vector line12 = point1 - point2;
const vector line23 = point2 - point3;
if
(
......@@ -87,141 +108,115 @@ void Foam::plane::calcPntAndVec
}
normal_ = line12 ^ line23;
scalar magUnitVector(mag(normal_));
if (magUnitVector < VSMALL)
{
FatalErrorInFunction
<< "Plane normal defined with zero length" << nl
<< "Bad points:" << point1 << ' ' << point2 << ' ' << point3
<< abort(FatalError);
}
normal_ /= magUnitVector;
makeUnitNormal(caller);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::plane::plane()
{}
Foam::plane::plane(const vector& normalVector)
:
normal_(normalVector),
point_(Zero)
origin_(Zero)
{
scalar magUnitVector(mag(normal_));
if (magUnitVector > VSMALL)
{
normal_ /= magUnitVector;
}
else
{
FatalErrorInFunction
<< "plane normal has zero length. base point:" << point_
<< abort(FatalError);
}
makeUnitNormal(FUNCTION_NAME);
}
Foam::plane::plane
(
const point& basePoint,
const point& originPoint,
const vector& normalVector,
const bool normalise
const bool doNormalise
)
:
normal_(normalVector),
point_(basePoint)
origin_(originPoint)
{
scalar magSqrNormalVector(magSqr(normal_));
makeUnitNormal(FUNCTION_NAME, doNormalise);
}
if (magSqrNormalVector < VSMALL)
{
FatalErrorInFunction
<< "plane normal has zero length. base point:" << point_
<< abort(FatalError);
}
if (normalise)
{
normal_ /= Foam::sqrt(magSqrNormalVector);
}
Foam::plane::plane(const scalarList& coeffs)
{
calcFromCoeffs
(
coeffs[0],
coeffs[1],
coeffs[2],
coeffs[3],
FUNCTION_NAME
);