Commit 557d9163 authored by graham's avatar graham
Browse files

ENH: face::nearestPointClassify, similar to triangle.

parent 6cdbd0ad
......@@ -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
(
......
......@@ -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();
......
......@@ -360,8 +360,8 @@ void Foam::treeDataTriSurface::findNearest
// scalar b(E0& E1);
// scalar c(E1& E1);
// // Get nearest point in s,t coordinates (s is along E0, t is along
// // E1)
// // Get nearest point in s,t coordinates (s is along E0, t
// // is along E1)
// scalar s;
// scalar t;
......
......@@ -2220,40 +2220,40 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
if (nearType == triPointRef::NONE)
{
vector sampleNearestVec = (sample - nearestPoint);
// Nearest to face interior. Use faceNormal to determine side
<<<<<<< HEAD
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.
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);
}
}
>>>>>>> 0bb6ebd... ENH: Making nearestPointClassify query for triangle.
// // 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)
{
......@@ -2273,26 +2273,26 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
label edgeI = surf.faceEdges()[nearestFaceI][nearLabel];
// 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
!= edge(f[nearLabel], f[f.fcIndex(nearLabel)])
)
{
FatalErrorIn("triSurfaceTools::surfaceSide")
<< "Edge:" << edgeI << " local vertices:" << e
<< " mesh vertices:" << meshEdge
<< " not at position " << nearLabel
<< " in face " << f
<< abort(FatalError);
}
}
// {
// // 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
// != edge(f[nearLabel], f[f.fcIndex(nearLabel)])
// )
// {
// FatalErrorIn("triSurfaceTools::surfaceSide")
// << "Edge:" << edgeI << " local vertices:" << e
// << " mesh vertices:" << meshEdge
// << " not at position " << nearLabel
// << " in face " << f
// << abort(FatalError);
// }
// }
return edgeSide(surf, sample, nearestPoint, edgeI);
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment