diff --git a/src/sampling/cuttingPlane/cuttingPlane.C b/src/sampling/cuttingPlane/cuttingPlane.C index f546c31dbe4ddfd5d5aed6d4998b392111c8ded5..472d30e7827f59297efc1f1442bc72c2020fd573 100644 --- a/src/sampling/cuttingPlane/cuttingPlane.C +++ b/src/sampling/cuttingPlane/cuttingPlane.C @@ -29,6 +29,15 @@ License #include "linePointRef.H" #include "meshTools.H" +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +// set values for what is close to zero and what is considered to +// be positive (and not just rounding noise) +//! @cond localScope +const Foam::scalar zeroish = Foam::SMALL; +const Foam::scalar positive = Foam::SMALL * 1E3; +//! @endcond localScope + // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // Find cut cells @@ -71,8 +80,8 @@ void Foam::cuttingPlane::calcCutCells if ( - (dotProducts[e[0]] < 0 && dotProducts[e[1]] > 0) - || (dotProducts[e[1]] < 0 && dotProducts[e[0]] > 0) + (dotProducts[e[0]] < zeroish && dotProducts[e[1]] > positive) + || (dotProducts[e[1]] < zeroish && dotProducts[e[0]] > positive) ) { nCutEdges++; @@ -116,8 +125,8 @@ Foam::labelList Foam::cuttingPlane::intersectEdges if ( - (dotProducts[e[0]] < 0 && dotProducts[e[1]] > 0) - || (dotProducts[e[1]] < 0 && dotProducts[e[0]] > 0) + (dotProducts[e[0]] < zeroish && dotProducts[e[1]] > positive) + || (dotProducts[e[1]] < zeroish && dotProducts[e[0]] > positive) ) { // Edge is cut. @@ -126,7 +135,19 @@ Foam::labelList Foam::cuttingPlane::intersectEdges scalar alpha = lineIntersect(linePointRef(p0, p1)); - dynCuttingPoints.append((1-alpha)*p0 + alpha*p1); + if (alpha < zeroish) + { + dynCuttingPoints.append(p0); + } + else if (alpha > 1.0) + { + dynCuttingPoints.append(p1); + } + else + { + dynCuttingPoints.append((1-alpha)*p0 + alpha*p1); + } + edgePoint[edgeI] = dynCuttingPoints.size() - 1; } } diff --git a/src/sampling/cuttingPlane/cuttingPlane.H b/src/sampling/cuttingPlane/cuttingPlane.H index 32dec6908d9c1248d1d3faedb4e1944b3a497221..3e33be9da769c887a842b0a7c7b61084cc0bb671 100644 --- a/src/sampling/cuttingPlane/cuttingPlane.H +++ b/src/sampling/cuttingPlane/cuttingPlane.H @@ -30,6 +30,10 @@ Description No attempt at resolving degenerate cases. +Note + When the cutting plane coincides with a mesh face, the cell edge on the + positive side of the plane is taken. + SourceFiles cuttingPlane.C