diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index d9548b67b10db88609ca0f132f6652b569fff356..94bc86f2efeacf8cfc767223e141b8b27c02f058 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -296,6 +296,7 @@ primitiveShapes = meshes/primitiveShapes
$(primitiveShapes)/line/line.C
$(primitiveShapes)/plane/plane.C
$(primitiveShapes)/triangle/intersection.C
+$(primitiveShapes)/objectHit/pointIndexHitIOList.C
meshShapes = meshes/meshShapes
$(meshShapes)/edge/edge.C
diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.H b/src/OpenFOAM/meshes/meshShapes/face/face.H
index 81ce90d8a9d537d2021ab3393bfb50235094f93c..dd5fc4a832b4265e0cd9aa0dd5842c224c13ba22 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/face.H
+++ b/src/OpenFOAM/meshes/meshShapes/face/face.H
@@ -126,6 +126,14 @@ class face
public:
+ //- Return types for classify
+ enum proxType
+ {
+ NONE,
+ POINT, // Close to point
+ EDGE // Close to edge
+ };
+
// Static data members
static const char* const typeName;
@@ -249,6 +257,20 @@ public:
const pointField& meshPoints
) const;
+ //- Return nearest point to face and classify it:
+ // + near point (nearType=POINT, nearLabel=0, 1, 2)
+ // + near edge (nearType=EDGE, nearLabel=0, 1, 2)
+ // Note: edges are counted from starting vertex so
+ // e.g. edge n is from f[n] to f[0], where the face has n + 1
+ // points
+ pointHit nearestPointClassify
+ (
+ const point& p,
+ const pointField& meshPoints,
+ label& nearType,
+ label& nearLabel
+ ) const;
+
//- Return contact sphere diameter
scalar contactSphereDiameter
(
diff --git a/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C b/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C
index d4184e766aa6a3dd225942da99c4d1e964330e83..f9f02eb61293307a29975ff878f75ffaef138ec0 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C
+++ b/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C
@@ -181,6 +181,22 @@ Foam::pointHit Foam::face::nearestPoint
const point& p,
const pointField& meshPoints
) const
+{
+ // Dummy labels
+ label nearType = -1;
+ label nearLabel = -1;
+
+ return nearestPointClassify(p, meshPoints, nearType, nearLabel);
+}
+
+
+Foam::pointHit Foam::face::nearestPointClassify
+(
+ const point& p,
+ const pointField& meshPoints,
+ label& nearType,
+ label& nearLabel
+) const
{
const face& f = *this;
point ctr = centre(meshPoints);
@@ -188,6 +204,9 @@ Foam::pointHit Foam::face::nearestPoint
// Initialize to miss, distance=GREAT
pointHit nearest(p);
+ nearType = -1;
+ nearLabel = -1;
+
label nPoints = f.size();
point nextPoint = ctr;
@@ -196,8 +215,10 @@ Foam::pointHit Foam::face::nearestPoint
{
nextPoint = meshPoints[f[fcIndex(pI)]];
+ label tmpNearType = -1;
+ label tmpNearLabel = -1;
+
// Note: for best accuracy, centre point always comes last
- //
triPointRef tri
(
meshPoints[f[pI]],
@@ -205,12 +226,42 @@ Foam::pointHit Foam::face::nearestPoint
ctr
);
- pointHit curHit = tri.nearestPoint(p);
+ pointHit curHit = tri.nearestPointClassify
+ (
+ p,
+ tmpNearType,
+ tmpNearLabel
+ );
if (Foam::mag(curHit.distance()) < Foam::mag(nearest.distance()))
{
nearest.setDistance(curHit.distance());
+ // Assume at first that the near type is NONE on the
+ // triangle (i.e. on the face of the triangle) then it is
+ // therefore also for the face.
+
+ nearType = NONE;
+
+ if (tmpNearType == triPointRef::EDGE && tmpNearLabel == 0)
+ {
+ // If the triangle edge label is 0, then this is also
+ // an edge of the face, if not, it is on the face
+
+ nearType = EDGE;
+
+ nearLabel = pI;
+ }
+ else if (tmpNearType == triPointRef::POINT && tmpNearLabel < 2)
+ {
+ // If the triangle point label is 0 or 1, then this is
+ // also a point of the face, if not, it is on the face
+
+ nearType = POINT;
+
+ nearLabel = pI + tmpNearLabel;
+ }
+
if (curHit.hit())
{
nearest.setHit();
diff --git a/src/meshTools/octree/PointIndexHit.H b/src/OpenFOAM/meshes/primitiveShapes/objectHit/PointIndexHit.H
similarity index 100%
rename from src/meshTools/octree/PointIndexHit.H
rename to src/OpenFOAM/meshes/primitiveShapes/objectHit/PointIndexHit.H
diff --git a/src/meshTools/octree/pointHitSort.H b/src/OpenFOAM/meshes/primitiveShapes/objectHit/pointHitSort.H
similarity index 100%
rename from src/meshTools/octree/pointHitSort.H
rename to src/OpenFOAM/meshes/primitiveShapes/objectHit/pointHitSort.H
diff --git a/src/meshTools/octree/pointIndexHit.H b/src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHit.H
similarity index 100%
rename from src/meshTools/octree/pointIndexHit.H
rename to src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHit.H
diff --git a/src/meshTools/octree/pointIndexHitIOList.C b/src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHitIOList.C
similarity index 100%
rename from src/meshTools/octree/pointIndexHitIOList.C
rename to src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHitIOList.C
diff --git a/src/meshTools/octree/pointIndexHitIOList.H b/src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHitIOList.H
similarity index 100%
rename from src/meshTools/octree/pointIndexHitIOList.H
rename to src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHitIOList.H
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
index 64d9318306fa0322ac3626fccf89bc5f46cc8be6..abb5d7d8e426ea94438536986ff92d9c731e3be0 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
@@ -79,21 +79,6 @@ class triangle
PointRef a_, b_, c_;
- // Private Member Functions
-
- //- Fast distance to triangle calculation. From
- // "Distance Between Point and Trangle in 3D"
- // David Eberly, Magic Software Inc. Aug. 2002.
- // Works on function Q giving distance to point and tries to
- // minimize this.
- static pointHit nearestPoint
- (
- const Point& baseVertex,
- const vector& E0,
- const vector& E1,
- const point& P
- );
-
public:
@@ -202,24 +187,27 @@ public:
const scalar tol = 0.0
) const;
- //- Return nearest point to p on triangle
- inline pointHit nearestPoint
+ //- Find the nearest point to p on the triangle and classify it:
+ // + near point (nearType=POINT, nearLabel=0, 1, 2)
+ // + near edge (nearType=EDGE, nearLabel=0, 1, 2)
+ // Note: edges are counted from starting
+ // vertex so e.g. edge 2 is from f[2] to f[0]
+ pointHit nearestPointClassify
(
- const point& p
+ const point& p,
+ label& nearType,
+ label& nearLabel
) const;
- //- Classify point in triangle plane w.r.t. triangle edges.
- // - inside (true returned)/outside (false returned)
- // - near point (nearType=POINT, nearLabel=0, 1, 2)
- // - near edge (nearType=EDGE, nearLabel=0, 1, 2)
- // Note: edges are counted from starting
- // vertex so e.g. edge 2 is from f[2] to f[0]
- // tol is fraction to account for truncation error. Is only used
- // when comparing normalized (0..1) numbers.
+ //- Return nearest point to p on triangle
+ inline pointHit nearestPoint(const point& p) const;
+
+ //- Classify nearest point to p in triangle plane
+ // w.r.t. triangle edges and points. Returns inside
+ // (true)/outside (false).
bool classify
(
const point& p,
- const scalar tol,
label& nearType,
label& nearLabel
) const;
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
index 10267a923f720824e5f50a7606f31ae62089fb07..7840d8ee7b4762c7e69c08b5238cb536a486a5a4 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
@@ -32,158 +32,6 @@ License
namespace Foam
{
-// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-
-template<class Point, class PointRef>
-pointHit triangle<Point, PointRef>::nearestPoint
-(
- const Point& baseVertex,
- const vector& E0,
- const vector& E1,
- const point& P
-)
-{
- // Distance vector
- const vector D(baseVertex - P);
-
- // Some geometrical factors
- const scalar a = E0 & E0;
- const scalar b = E0 & E1;
- const scalar c = E1 & E1;
-
- // Precalculate distance factors
- const scalar d = E0 & D;
- const scalar e = E1 & D;
- const scalar f = D & D;
-
- // Do classification
- const scalar det = a*c - b*b;
- scalar s = b*e - c*d;
- scalar t = b*d - a*e;
-
- bool inside = false;
-
- if (s+t < det)
- {
- if (s < 0)
- {
- if (t < 0)
- {
- // Region 4
- if (e > 0)
- {
- // min on edge t = 0
- t = 0;
- s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
- }
- else
- {
- // min on edge s=0
- s = 0;
- t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
- }
- }
- else
- {
- // Region 3. Min on edge s = 0
- s = 0;
- t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
- }
- }
- else if (t < 0)
- {
- // Region 5
- t = 0;
- s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
- }
- else
- {
- // Region 0
- const scalar invDet = 1/det;
- s *= invDet;
- t *= invDet;
-
- inside = true;
- }
- }
- else
- {
- if (s < 0)
- {
- // Region 2
- const scalar tmp0 = b + d;
- const scalar tmp1 = c + e;
- if (tmp1 > tmp0)
- {
- // min on edge s+t=1
- const scalar numer = tmp1 - tmp0;
- const scalar denom = a-2*b+c;
- s = (numer >= denom ? 1 : numer/denom);
- t = 1 - s;
- }
- else
- {
- // min on edge s=0
- s = 0;
- t = (tmp1 <= 0 ? 1 : (e >= 0 ? 0 : - e/c));
- }
- }
- else if (t < 0)
- {
- // Region 6
- const scalar tmp0 = b + d;
- const scalar tmp1 = c + e;
- if (tmp1 > tmp0)
- {
- // min on edge s+t=1
- const scalar numer = tmp1 - tmp0;
- const scalar denom = a-2*b+c;
- s = (numer >= denom ? 1 : numer/denom);
- t = 1 - s;
- }
- else
- {
- // min on edge t=0
- t = 0;
- s = (tmp1 <= 0 ? 1 : (d >= 0 ? 0 : - d/a));
- }
- }
- else
- {
- // Region 1
- const scalar numer = c+e-(b+d);
- if (numer <= 0)
- {
- s = 0;
- }
- else
- {
- const scalar denom = a-2*b+c;
- s = (numer >= denom ? 1 : numer/denom);
- }
- }
-
- t = 1 - s;
- }
-
- // Calculate distance.
- // Note: Foam::mag used since truncation error causes negative distances
- // with points very close to one of the triangle vertices.
- // (Up to -2.77556e-14 seen). Could use +SMALL but that not large enough.
-
- return pointHit
- (
- inside,
- baseVertex + s*E0 + t*E1,
- Foam::sqrt
- (
- Foam::mag(a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f)
- ),
- !inside
- );
-}
-
-
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Point, class PointRef>
@@ -247,7 +95,7 @@ inline Point triangle<Point, PointRef>::centre() const
template<class Point, class PointRef>
inline scalar triangle<Point, PointRef>::mag() const
{
- return ::Foam::mag(normal());
+ return Foam::mag(normal());
}
@@ -536,7 +384,7 @@ inline pointHit triangle<Point, PointRef>::ray
inter.setMiss(eligible);
// The miss point is the nearest point on the triangle
- inter.setPoint(nearestPoint(a_, E0, E1, p).rawPoint());
+ inter.setPoint(nearestPoint(p).rawPoint());
// The distance to the miss is the distance between the
// original point and plane of intersection
@@ -634,177 +482,143 @@ inline pointHit triangle<Point, PointRef>::intersection
template<class Point, class PointRef>
-inline pointHit triangle<Point, PointRef>::nearestPoint
-(
- const point& p
-) const
-{
- // Express triangle in terms of baseVertex (point a_) and
- // two edge vectors
- vector E0 = b_ - a_;
- vector E1 = c_ - a_;
-
- return nearestPoint(a_, E0, E1, p);
-}
-
-
-template<class Point, class PointRef>
-inline bool triangle<Point, PointRef>::classify
+pointHit triangle<Point, PointRef>::nearestPointClassify
(
const point& p,
- const scalar tol,
label& nearType,
label& nearLabel
) const
{
- const vector E0 = b_ - a_;
- const vector E1 = c_ - a_;
- const vector n = 0.5*(E0 ^ E1);
-
- // Get largest component of normal
- scalar magX = Foam::mag(n.x());
- scalar magY = Foam::mag(n.y());
- scalar magZ = Foam::mag(n.z());
-
- label i0 = -1;
- if ((magX >= magY) && (magX >= magZ))
- {
- i0 = 0;
- }
- else if ((magY >= magX) && (magY >= magZ))
- {
- i0 = 1;
- }
- else
- {
- i0 = 2;
- }
+ // Adapted from:
+ // Real-time collision detection, Christer Ericson, 2005, 136-142
- // Get other components
- label i1 = (i0 + 1) % 3;
- label i2 = (i1 + 1) % 3;
+ // Check if P in vertex region outside A
+ vector ab = b_ - a_;
+ vector ac = c_ - a_;
+ vector ap = p - a_;
+ scalar d1 = ab & ap;
+ scalar d2 = ac & ap;
- scalar u1 = E0[i1];
- scalar v1 = E0[i2];
+ if (d1 <= 0.0 && d2 <= 0.0)
+ {
+ // barycentric coordinates (1, 0, 0)
- scalar u2 = E1[i1];
- scalar v2 = E1[i2];
+ nearType = POINT;
+ nearLabel = 0;
+ return pointHit(false, a_, Foam::mag(a_ - p), true);
+ }
- scalar det = v2*u1 - u2*v1;
+ // Check if P in vertex region outside B
+ vector bp = p - b_;
+ scalar d3 = ab & bp;
+ scalar d4 = ac & bp;
- scalar u0 = p[i1] - a_[i1];
- scalar v0 = p[i2] - a_[i2];
+ if (d3 >= 0.0 && d4 <= d3)
+ {
+ // barycentric coordinates (0, 1, 0)
- scalar alpha = 0;
- scalar beta = 0;
+ nearType = POINT;
+ nearLabel = 1;
+ return pointHit(false, b_, Foam::mag(b_ - p), true);
+ }
- bool hit = false;
+ // Check if P in edge region of AB, if so return projection of P onto AB
+ scalar vc = d1*d4 - d3*d2;
- if (Foam::mag(u1) < ROOTVSMALL)
+ if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0)
{
- beta = u0/u2;
-
- alpha = (v0 - beta*v2)/v1;
+ // barycentric coordinates (1-v, v, 0)
+ scalar v = d1/(d1 - d3);
- hit = ((alpha >= 0) && ((alpha + beta) <= 1));
+ point nearPt = a_ + v*ab;
+ nearType = EDGE;
+ nearLabel = 0;
+ return pointHit(false, nearPt, Foam::mag(nearPt - p), true);
}
- else
- {
- beta = (v0*u1 - u0*v1)/det;
- alpha = (u0 - beta*u2)/u1;
+ // Check if P in vertex region outside C
+ vector cp = p - c_;
+ scalar d5 = ab & cp;
+ scalar d6 = ac & cp;
- hit = ((alpha >= 0) && ((alpha + beta) <= 1));
- }
+ if (d6 >= 0.0 && d5 <= d6)
+ {
+ // barycentric coordinates (0, 0, 1)
- //
- // Now alpha, beta are the coordinates in the triangle local coordinate
- // system E0, E1
- //
+ nearType = POINT;
+ nearLabel = 2;
+ return pointHit(false, c_, Foam::mag(c_ - p), true);
+ }
- //Pout<< "alpha:" << alpha << endl;
- //Pout<< "beta:" << beta << endl;
- //Pout<< "hit:" << hit << endl;
- //Pout<< "tol:" << tol << endl;
+ // Check if P in edge region of AC, if so return projection of P onto AC
+ scalar vb = d5*d2 - d1*d6;
- if (hit)
+ if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0)
{
- // alpha,beta might get negative due to precision errors
- alpha = max(0.0, min(1.0, alpha));
- beta = max(0.0, min(1.0, beta));
+ // barycentric coordinates (1-w, 0, w)
+ scalar w = d2/(d2 - d6);
+
+ point nearPt = a_ + w*ac;
+ nearType = EDGE;
+ nearLabel = 2;
+ return pointHit(false, nearPt, Foam::mag(nearPt - p), true);
}
- nearType = NONE;
- nearLabel = -1;
+ // Check if P in edge region of BC, if so return projection of P onto BC
+ scalar va = d3*d6 - d5*d4;
- if (Foam::mag(alpha+beta-1) <= tol)
+ if (va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0)
{
- // On edge between vert 1-2 (=edge 1)
+ // barycentric coordinates (0, 1-w, w)
+ scalar w = (d4 - d3)/((d4 - d3) + (d5 - d6));
- if (Foam::mag(alpha) <= tol)
- {
- nearType = POINT;
- nearLabel = 2;
- }
- else if (Foam::mag(beta) <= tol)
- {
- nearType = POINT;
- nearLabel = 1;
- }
- else if ((alpha >= 0) && (alpha <= 1) && (beta >= 0) && (beta <= 1))
- {
- nearType = EDGE;
- nearLabel = 1;
- }
+ point nearPt = b_ + w*(c_ - b_);
+ nearType = EDGE;
+ nearLabel = 1;
+ return pointHit(false, nearPt, Foam::mag(nearPt - p), true);
}
- else if (Foam::mag(alpha) <= tol)
- {
- // On edge between vert 2-0 (=edge 2)
- if (Foam::mag(beta) <= tol)
- {
- nearType = POINT;
- nearLabel = 0;
- }
- else if (Foam::mag(beta-1) <= tol)
- {
- nearType = POINT;
- nearLabel = 2;
- }
- else if ((beta >= 0) && (beta <= 1))
- {
- nearType = EDGE;
- nearLabel = 2;
- }
- }
- else if (Foam::mag(beta) <= tol)
- {
- // On edge between vert 0-1 (= edge 0)
+ // P inside face region. Compute Q through its barycentric
+ // coordinates (u, v, w)
+ scalar denom = 1.0/(va + vb + vc);
+ scalar v = vb * denom;
+ scalar w = vc * denom;
- if (Foam::mag(alpha) <= tol)
- {
- nearType = POINT;
- nearLabel = 0;
- }
- else if (Foam::mag(alpha-1) <= tol)
- {
- nearType = POINT;
- nearLabel = 1;
- }
- else if ((alpha >= 0) && (alpha <= 1))
- {
- nearType = EDGE;
- nearLabel = 0;
- }
- }
+ // = u*a + v*b + w*c, u = va*denom = 1.0 - v - w
- return hit;
+ point nearPt = a_ + ab*v + ac*w;
+ nearType = NONE,
+ nearLabel = -1;
+ return pointHit(true, nearPt, Foam::mag(nearPt - p), false);
}
+template<class Point, class PointRef>
+inline pointHit triangle<Point, PointRef>::nearestPoint
+(
+ const point& p
+) const
+{
+ // Dummy labels
+ label nearType = -1;
+ label nearLabel = -1;
+
+ return nearestPointClassify(p, nearType, nearLabel);
+}
+template<class Point, class PointRef>
+inline bool triangle<Point, PointRef>::classify
+(
+ const point& p,
+ label& nearType,
+ label& nearLabel
+) const
+{
+ return nearestPointClassify(p, nearType, nearLabel).hit();
+}
+
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index 14e561fcbab7221adc85de4628a2beb592841680..2c9f81195d082b75e3b08a62ede15f2a8ec37d85 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -44,7 +44,6 @@ octree/octreeDataFace.C
octree/treeBoundBox.C
octree/treeNodeName.C
octree/treeLeafName.C
-octree/pointIndexHitIOList.C
indexedOctree/indexedOctreeName.C
indexedOctree/treeDataCell.C
diff --git a/src/meshTools/indexedOctree/treeDataTriSurface.C b/src/meshTools/indexedOctree/treeDataTriSurface.C
index e0f6d1a9df32cd37591245158e4fbce0620d0878..f380167174481939ae61a056bae233485239d815 100644
--- a/src/meshTools/indexedOctree/treeDataTriSurface.C
+++ b/src/meshTools/indexedOctree/treeDataTriSurface.C
@@ -35,145 +35,145 @@ defineTypeNameAndDebug(Foam::treeDataTriSurface, 0);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-// Fast distance to triangle calculation. From
-// "Distance Between Point and Trangle in 3D"
-// David Eberly, Magic Software Inc. Aug. 2003.
-// Works on function Q giving distance to point and tries to minimize this.
-Foam::scalar Foam::treeDataTriSurface::nearestCoords
-(
- const point& base,
- const point& E0,
- const point& E1,
- const scalar a,
- const scalar b,
- const scalar c,
- const point& P,
- scalar& s,
- scalar& t
-)
-{
- // distance vector
- const vector D(base - P);
-
- // Precalculate distance factors.
- const scalar d = E0 & D;
- const scalar e = E1 & D;
-
- // Do classification
- const scalar det = a*c - b*b;
-
- s = b*e - c*d;
- t = b*d - a*e;
-
- if (s+t < det)
- {
- if (s < 0)
- {
- if (t < 0)
- {
- //region 4
- if (e > 0)
- {
- //min on edge t = 0
- t = 0;
- s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
- }
- else
- {
- //min on edge s=0
- s = 0;
- t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
- }
- }
- else
- {
- //region 3. Min on edge s = 0
- s = 0;
- t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
- }
- }
- else if (t < 0)
- {
- //region 5
- t = 0;
- s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
- }
- else
- {
- //region 0
- const scalar invDet = 1/det;
- s *= invDet;
- t *= invDet;
- }
- }
- else
- {
- if (s < 0)
- {
- //region 2
- const scalar tmp0 = b + d;
- const scalar tmp1 = c + e;
- if (tmp1 > tmp0)
- {
- //min on edge s+t=1
- const scalar numer = tmp1 - tmp0;
- const scalar denom = a-2*b+c;
- s = (numer >= denom ? 1 : numer/denom);
- t = 1 - s;
- }
- else
- {
- //min on edge s=0
- s = 0;
- t = (tmp1 <= 0 ? 1 : (e >= 0 ? 0 : - e/c));
- }
- }
- else if (t < 0)
- {
- //region 6
- const scalar tmp0 = b + d;
- const scalar tmp1 = c + e;
- if (tmp1 > tmp0)
- {
- //min on edge s+t=1
- const scalar numer = tmp1 - tmp0;
- const scalar denom = a-2*b+c;
- s = (numer >= denom ? 1 : numer/denom);
- t = 1 - s;
- }
- else
- {
- //min on edge t=0
- t = 0;
- s = (tmp1 <= 0 ? 1 : (d >= 0 ? 0 : - d/a));
- }
- }
- else
- {
- //region 1
- const scalar numer = c+e-(b+d);
- if (numer <= 0)
- {
- s = 0;
- }
- else
- {
- const scalar denom = a-2*b+c;
- s = (numer >= denom ? 1 : numer/denom);
- }
- }
- t = 1 - s;
- }
-
-
- // Calculate distance.
- // Note: abs should not be needed but truncation error causes problems
- // with points very close to one of the triangle vertices.
- // (seen up to -9e-15). Alternatively add some small value.
-
- const scalar f = D & D;
- return Foam::mag(a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f);
-}
+// // Fast distance to triangle calculation. From
+// // "Distance Between Point and Triangle in 3D"
+// // David Eberly, Magic Software Inc. Aug. 2003.
+// // Works on function Q giving distance to point and tries to minimize this.
+// Foam::scalar Foam::treeDataTriSurface::nearestCoords
+// (
+// const point& base,
+// const point& E0,
+// const point& E1,
+// const scalar a,
+// const scalar b,
+// const scalar c,
+// const point& P,
+// scalar& s,
+// scalar& t
+// )
+// {
+// // distance vector
+// const vector D(base - P);
+
+// // Precalculate distance factors.
+// const scalar d = E0 & D;
+// const scalar e = E1 & D;
+
+// // Do classification
+// const scalar det = a*c - b*b;
+
+// s = b*e - c*d;
+// t = b*d - a*e;
+
+// if (s+t < det)
+// {
+// if (s < 0)
+// {
+// if (t < 0)
+// {
+// //region 4
+// if (e > 0)
+// {
+// //min on edge t = 0
+// t = 0;
+// s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
+// }
+// else
+// {
+// //min on edge s=0
+// s = 0;
+// t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
+// }
+// }
+// else
+// {
+// //region 3. Min on edge s = 0
+// s = 0;
+// t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
+// }
+// }
+// else if (t < 0)
+// {
+// //region 5
+// t = 0;
+// s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
+// }
+// else
+// {
+// //region 0
+// const scalar invDet = 1/det;
+// s *= invDet;
+// t *= invDet;
+// }
+// }
+// else
+// {
+// if (s < 0)
+// {
+// //region 2
+// const scalar tmp0 = b + d;
+// const scalar tmp1 = c + e;
+// if (tmp1 > tmp0)
+// {
+// //min on edge s+t=1
+// const scalar numer = tmp1 - tmp0;
+// const scalar denom = a-2*b+c;
+// s = (numer >= denom ? 1 : numer/denom);
+// t = 1 - s;
+// }
+// else
+// {
+// //min on edge s=0
+// s = 0;
+// t = (tmp1 <= 0 ? 1 : (e >= 0 ? 0 : - e/c));
+// }
+// }
+// else if (t < 0)
+// {
+// //region 6
+// const scalar tmp0 = b + d;
+// const scalar tmp1 = c + e;
+// if (tmp1 > tmp0)
+// {
+// //min on edge s+t=1
+// const scalar numer = tmp1 - tmp0;
+// const scalar denom = a-2*b+c;
+// s = (numer >= denom ? 1 : numer/denom);
+// t = 1 - s;
+// }
+// else
+// {
+// //min on edge t=0
+// t = 0;
+// s = (tmp1 <= 0 ? 1 : (d >= 0 ? 0 : - d/a));
+// }
+// }
+// else
+// {
+// //region 1
+// const scalar numer = c+e-(b+d);
+// if (numer <= 0)
+// {
+// s = 0;
+// }
+// else
+// {
+// const scalar denom = a-2*b+c;
+// s = (numer >= denom ? 1 : numer/denom);
+// }
+// }
+// t = 1 - s;
+// }
+
+
+// // Calculate distance.
+// // Note: abs should not be needed but truncation error causes problems
+// // with points very close to one of the triangle vertices.
+// // (seen up to -9e-15). Alternatively add some small value.
+
+// const scalar f = D & D;
+// return Foam::mag(a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f);
+// }
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@@ -234,9 +234,7 @@ Foam::label Foam::treeDataTriSurface::getVolumeType
(
surface_,
sample,
- pHit.index(),
- pHit.hitPoint(),
- indexedOctree<treeDataTriSurface>::perturbTol()
+ pHit.index()
);
if (t == triSurfaceTools::UNKNOWN)
@@ -353,39 +351,43 @@ void Foam::treeDataTriSurface::findNearest
// )
//)
{
- // Get spanning vectors of triangle
- vector base(p1);
- vector E0(p0 - p1);
- vector E1(p2 - p1);
-
- scalar a(E0& E0);
- scalar b(E0& E1);
- scalar c(E1& E1);
-
- // Get nearest point in s,t coordinates (s is along E0, t is along
- // E1)
- scalar s;
- scalar t;
-
- scalar distSqr = nearestCoords
- (
- base,
- E0,
- E1,
- a,
- b,
- c,
- sample,
-
- s,
- t
- );
+ // // Get spanning vectors of triangle
+ // vector base(p1);
+ // vector E0(p0 - p1);
+ // vector E1(p2 - p1);
+
+ // scalar a(E0& E0);
+ // scalar b(E0& E1);
+ // scalar c(E1& E1);
+
+ // // Get nearest point in s,t coordinates (s is along E0, t
+ // // is along E1)
+ // scalar s;
+ // scalar t;
+
+ // scalar distSqr = nearestCoords
+ // (
+ // base,
+ // E0,
+ // E1,
+ // a,
+ // b,
+ // c,
+ // sample,
+
+ // s,
+ // t
+ // );
+
+ pointHit pHit = triPointRef(p0, p1, p2).nearestPoint(sample);
+
+ scalar distSqr = sqr(pHit.distance());
if (distSqr < nearestDistSqr)
{
nearestDistSqr = distSqr;
minIndex = index;
- nearestPoint = base + s*E0 + t*E1;
+ nearestPoint = pHit.rawPoint();
}
}
}
diff --git a/src/meshTools/indexedOctree/treeDataTriSurface.H b/src/meshTools/indexedOctree/treeDataTriSurface.H
index f602783afefe17913cd2484dd631f587245fd80a..dc959f78ff69ec71bb0f201d5fe9366f9afcc5a1 100644
--- a/src/meshTools/indexedOctree/treeDataTriSurface.H
+++ b/src/meshTools/indexedOctree/treeDataTriSurface.H
@@ -60,20 +60,20 @@ class treeDataTriSurface
// Private Member Functions
- //- fast triangle nearest point calculation. Returns point in E0, E1
- // coordinate system: base + s*E0 + t*E1
- static scalar nearestCoords
- (
- const point& base,
- const point& E0,
- const point& E1,
- const scalar a,
- const scalar b,
- const scalar c,
- const point& P,
- scalar& s,
- scalar& t
- );
+ // //- fast triangle nearest point calculation. Returns point in E0, E1
+ // // coordinate system: base + s*E0 + t*E1
+ // static scalar nearestCoords
+ // (
+ // const point& base,
+ // const point& E0,
+ // const point& E1,
+ // const scalar a,
+ // const scalar b,
+ // const scalar c,
+ // const point& P,
+ // scalar& s,
+ // scalar& t
+ // );
public:
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C
index 110f9dcd66161da1e9a66a3be1c711f718106bdb..c8955f747091e3b79bd7bbb8d0a30d8be4753b6f 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C
@@ -430,10 +430,7 @@ bool Foam::edgeIntersections::offsetPerturb
point ctr = tri.centre();
- // Get measure for tolerance.
- scalar tolDim = 0.001*mag(tri.a() - ctr);
-
- tri.classify(pHit.hitPoint(), tolDim, nearType, nearLabel);
+ tri.classify(pHit.hitPoint(), nearType, nearLabel);
if (nearType == triPointRef::POINT || nearType == triPointRef::EDGE)
{
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
index 6602fbcf430385f0018ae9368381c946571cdd1f..58c206aa75069667f64c5f55a738f452599084b8 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
@@ -315,7 +315,7 @@ void Foam::surfaceIntersection::classifyHit
surf2Pts[f2[0]],
surf2Pts[f2[1]],
surf2Pts[f2[2]]
- ).classify(pHit.hitPoint(), tolDim, nearType, nearLabel);
+ ).classify(pHit.hitPoint(), nearType, nearLabel);
// Classify points on edge of surface1
label edgeEnd =
diff --git a/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C b/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C
index c80991819ea2244fc4c1a1159f5a2e1836db30a8..55e60a3111cd5bd5c698a592edd70b5753ae603e 100644
--- a/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C
+++ b/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C
@@ -43,7 +43,7 @@ Foam::scalar Foam::octreeDataTriSurface::tol(1E-6);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Fast distance to triangle calculation. From
-// "Distance Between Point and Trangle in 3D"
+// "Distance Between Point and Triangle in 3D"
// David Eberly, Magic Software Inc. Aug. 2003.
// Works on function Q giving distance to point and tries to minimize this.
void Foam::octreeDataTriSurface::nearestCoords
diff --git a/src/meshTools/triSurface/orientedSurface/orientedSurface.C b/src/meshTools/triSurface/orientedSurface/orientedSurface.C
index 6462df4cb31ec2e08f197c5ddf1de6b6284597c3..38a7afd4b8ed5ad7524dcdc223a303695d672c2f 100644
--- a/src/meshTools/triSurface/orientedSurface/orientedSurface.C
+++ b/src/meshTools/triSurface/orientedSurface/orientedSurface.C
@@ -234,9 +234,7 @@ void Foam::orientedSurface::propagateOrientation
(
s,
samplePoint,
- nearestFaceI,
- nearestPt,
- 10*SMALL
+ nearestFaceI
);
if (side == triSurfaceTools::UNKNOWN)
diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
index 5a0cc0362340e7a36b5149a3ec721d09882653e2..af30b3bcc8954ad5e86916ff97d9d62e04c113cf 100644
--- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
+++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
@@ -2121,12 +2121,13 @@ Foam::vector Foam::triSurfaceTools::surfaceNormal
label nearType;
label nearLabel;
+
triPointRef
(
points[f[0]],
points[f[1]],
points[f[2]]
- ).classify(nearestPt, 1E-6, nearType, nearLabel);
+ ).classify(nearestPt, nearType, nearLabel);
if (nearType == triPointRef::NONE)
{
@@ -2199,28 +2200,61 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
(
const triSurface& surf,
const point& sample,
- const label nearestFaceI, // nearest face
- const point& nearestPoint, // nearest point on nearest face
- const scalar tol
+ const label nearestFaceI
)
{
const labelledTri& f = surf[nearestFaceI];
const pointField& points = surf.points();
- // Find where point is on triangle. Note tolerance needed. Is relative
- // one (used in comparing normalized [0..1] triangle coordinates).
+ // Find where point is on triangle.
label nearType, nearLabel;
- triPointRef
+
+ pointHit pHit = triPointRef
(
points[f[0]],
points[f[1]],
points[f[2]]
- ).classify(nearestPoint, tol, nearType, nearLabel);
+ ).nearestPointClassify(sample, nearType, nearLabel);
+
+ const point& nearestPoint(pHit.rawPoint());
if (nearType == triPointRef::NONE)
{
+ vector sampleNearestVec = (sample - nearestPoint);
+
// Nearest to face interior. Use faceNormal to determine side
- scalar c = (sample - nearestPoint) & surf.faceNormals()[nearestFaceI];
+ scalar c = sampleNearestVec & surf.faceNormals()[nearestFaceI];
+
+ // // If the sample is essentially on the face, do not check for
+ // // it being perpendicular.
+
+ // scalar magSampleNearestVec = mag(sampleNearestVec);
+
+ // if (magSampleNearestVec > SMALL)
+ // {
+ // c /= magSampleNearestVec*mag(surf.faceNormals()[nearestFaceI]);
+
+ // if (mag(c) < 0.99)
+ // {
+ // FatalErrorIn("triSurfaceTools::surfaceSide")
+ // << "nearestPoint identified as being on triangle face "
+ // << "but vector from nearestPoint to sample is not "
+ // << "perpendicular to the normal." << nl
+ // << "sample: " << sample << nl
+ // << "nearestPoint: " << nearestPoint << nl
+ // << "sample - nearestPoint: "
+ // << sample - nearestPoint << nl
+ // << "normal: " << surf.faceNormals()[nearestFaceI] << nl
+ // << "mag(sample - nearestPoint): "
+ // << mag(sample - nearestPoint) << nl
+ // << "normalised dot product: " << c << nl
+ // << "triangle vertices: " << nl
+ // << " " << points[f[0]] << nl
+ // << " " << points[f[1]] << nl
+ // << " " << points[f[2]] << nl
+ // << abort(FatalError);
+ // }
+ // }
if (c > 0)
{
@@ -2239,13 +2273,13 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
// Get the edge. Assume order of faceEdges same as face vertices.
label edgeI = surf.faceEdges()[nearestFaceI][nearLabel];
- //if (debug)
- //{
+ // if (debug)
+ // {
// // Check order of faceEdges same as face vertices.
// const edge& e = surf.edges()[edgeI];
// const labelList& meshPoints = surf.meshPoints();
// const edge meshEdge(meshPoints[e[0]], meshPoints[e[1]]);
- //
+
// if
// (
// meshEdge
@@ -2259,7 +2293,7 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
// << " in face " << f
// << abort(FatalError);
// }
- //}
+ // }
return edgeSide(surf, sample, nearestPoint, edgeI);
}
@@ -2717,7 +2751,14 @@ void Foam::triSurfaceTools::calcInterpolationWeights
triPointRef tri(f.tri(points));
- pointHit nearest = tri.nearestPoint(samplePt);
+ label nearType, nearLabel;
+
+ pointHit nearest = tri.nearestPointClassify
+ (
+ samplePt,
+ nearType,
+ nearLabel
+ );
if (nearest.hit())
{
@@ -2741,14 +2782,6 @@ void Foam::triSurfaceTools::calcInterpolationWeights
minDistance = nearest.distance();
// Outside triangle. Store nearest.
- label nearType, nearLabel;
- tri.classify
- (
- nearest.rawPoint(),
- 1E-6, // relative tolerance
- nearType,
- nearLabel
- );
if (nearType == triPointRef::POINT)
{
@@ -2779,12 +2812,12 @@ void Foam::triSurfaceTools::calcInterpolationWeights
max
(
0,
- mag(nearest.rawPoint()-p0)/mag(p1-p0)
+ mag(nearest.rawPoint() - p0)/mag(p1 - p0)
)
);
// Interpolate
- weights[0] = 1-s;
+ weights[0] = 1 - s;
weights[1] = s;
weights[2] = -GREAT;
@@ -2830,7 +2863,6 @@ Foam::surfaceLocation Foam::triSurfaceTools::classify
triPointRef(s[triI].tri(s.points())).classify
(
trianglePoint,
- 1E-6,
elemType,
index
);
diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H
index 62d66d031f02e5eb9e8723c8ef0f42b7d7bd0b87..744a12217d372adab1f4a6915d7c8fb3e7bfed29 100644
--- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H
+++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H
@@ -458,9 +458,7 @@ public:
(
const triSurface& surf,
const point& sample,
- const label nearestFaceI, // nearest face
- const point& nearestPt, // nearest point on nearest face
- const scalar tol // tolerance for nearness test.
+ const label nearestFaceI
);
// Triangulation of faces
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
index b777edbdb3c1a0d99713ed601ba0f664e9e313d3..4a4c9219d690dcb106b3b409f6994f535d5d65f1 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
@@ -89,20 +89,29 @@ void Foam::distanceSurface::createGeometry()
if (signed_)
{
- vectorField normal;
- surfPtr_().getNormal(nearest, normal);
+ List<searchableSurface::volumeType> volType;
- forAll(nearest, i)
+ surfPtr_().getVolumeType(cc, volType);
+
+ forAll(volType, i)
{
- vector d(cc[i]-nearest[i].hitPoint());
+ searchableSurface::volumeType vT = volType[i];
- if ((d&normal[i]) > 0)
+ if (vT == searchableSurface::OUTSIDE)
{
- fld[i] = Foam::mag(d);
+ fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
+ }
+ else if (vT == searchableSurface::INSIDE)
+ {
+ fld[i] = -Foam::mag(cc[i] - nearest[i].hitPoint());
}
else
{
- fld[i] = -Foam::mag(d);
+ FatalErrorIn
+ (
+ "void Foam::distanceSurface::createGeometry()"
+ ) << "getVolumeType failure, neither INSIDE or OUTSIDE"
+ << exit(FatalError);
}
}
}
@@ -132,20 +141,30 @@ void Foam::distanceSurface::createGeometry()
if (signed_)
{
- vectorField normal;
- surfPtr_().getNormal(nearest, normal);
+ List<searchableSurface::volumeType> volType;
- forAll(nearest, i)
+ surfPtr_().getVolumeType(cc, volType);
+
+ forAll(volType, i)
{
- vector d(cc[i]-nearest[i].hitPoint());
+ searchableSurface::volumeType vT = volType[i];
- if ((d&normal[i]) > 0)
+ if (vT == searchableSurface::OUTSIDE)
{
- fld[i] = Foam::mag(d);
+ fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
+ }
+ else if (vT == searchableSurface::INSIDE)
+ {
+ fld[i] = -Foam::mag(cc[i] - nearest[i].hitPoint());
}
else
{
- fld[i] = -Foam::mag(d);
+ FatalErrorIn
+ (
+ "void Foam::distanceSurface::createGeometry()"
+ ) << "getVolumeType failure, "
+ << "neither INSIDE or OUTSIDE"
+ << exit(FatalError);
}
}
}
@@ -179,20 +198,31 @@ void Foam::distanceSurface::createGeometry()
if (signed_)
{
- vectorField normal;
- surfPtr_().getNormal(nearest, normal);
+ List<searchableSurface::volumeType> volType;
- forAll(nearest, i)
+ surfPtr_().getVolumeType(pts, volType);
+
+ forAll(volType, i)
{
- vector d(pts[i]-nearest[i].hitPoint());
+ searchableSurface::volumeType vT = volType[i];
- if ((d&normal[i]) > 0)
+ if (vT == searchableSurface::OUTSIDE)
+ {
+ pointDistance_[i] =
+ Foam::mag(pts[i] - nearest[i].hitPoint());
+ }
+ else if (vT == searchableSurface::INSIDE)
{
- pointDistance_[i] = Foam::mag(d);
+ pointDistance_[i] =
+ -Foam::mag(pts[i] - nearest[i].hitPoint());
}
else
{
- pointDistance_[i] = -Foam::mag(d);
+ FatalErrorIn
+ (
+ "void Foam::distanceSurface::createGeometry()"
+ ) << "getVolumeType failure, neither INSIDE or OUTSIDE"
+ << exit(FatalError);
}
}
}