diff --git a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
index 024f7d27ecf4203c1684e3295bc7e20ad52692fc..f16dd2595c523183104a52a8330d62588e8f2758 100644
--- a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
+++ b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
@@ -2832,7 +2832,7 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
) const
{
label nearestShapeI = -1;
- point nearestPoint;
+ point nearestPoint = vector::zero;
if (nodes_.size())
{
@@ -2847,10 +2847,6 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
nearestPoint
);
}
- else
- {
- nearestPoint = vector::zero;
- }
return pointIndexHit(nearestShapeI != -1, nearestPoint, nearestShapeI);
}
diff --git a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
index e39eae5704a61cbc53438060ac1a8ca77154144e..465734ff879bd4aa2507b05fe3779cb1d4984ea9 100644
--- a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
+++ b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
@@ -33,11 +33,15 @@ License
#include "OFstream.H"
#include "IFstream.H"
#include "unitConversion.H"
+#include "EdgeMap.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::surfaceFeatures, 0);
+const Foam::scalar Foam::surfaceFeatures::parallelTolerance =
+ sin(degToRad(1.0));
+
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@@ -170,7 +174,10 @@ void Foam::surfaceFeatures::setFromStatus(const List<edgeStatus>& edgeStat)
//construct feature points where more than 2 feature edges meet
-void Foam::surfaceFeatures::calcFeatPoints(const List<edgeStatus>& edgeStat)
+void Foam::surfaceFeatures::calcFeatPoints
+(
+ const List<edgeStatus>& edgeStat
+)
{
DynamicList<label> featurePoints(surf_.nPoints()/1000);
@@ -200,6 +207,60 @@ void Foam::surfaceFeatures::calcFeatPoints(const List<edgeStatus>& edgeStat)
}
+void Foam::surfaceFeatures::classifyFeatureAngles
+(
+ const labelListList& edgeFaces,
+ List<edgeStatus>& edgeStat,
+ const scalar minCos
+) const
+{
+ const vectorField& faceNormals = surf_.faceNormals();
+ const pointField& points = surf_.points();
+
+ forAll(edgeFaces, edgeI)
+ {
+ const labelList& eFaces = edgeFaces[edgeI];
+
+ if (eFaces.size() != 2)
+ {
+ // Non-manifold. What to do here? Is region edge? external edge?
+ edgeStat[edgeI] = REGION;
+ }
+ else
+ {
+ label face0 = eFaces[0];
+ label face1 = eFaces[1];
+
+ if (surf_[face0].region() != surf_[face1].region())
+ {
+ edgeStat[edgeI] = REGION;
+ }
+ else
+ {
+ if ((faceNormals[face0] & faceNormals[face1]) < minCos)
+ {
+
+ // Check if convex or concave by looking at angle
+ // between face centres and normal
+ vector f0Tof1 =
+ surf_[face1].centre(points)
+ - surf_[face0].centre(points);
+
+ if ((f0Tof1 & faceNormals[face0]) > 0.0)
+ {
+ edgeStat[edgeI] = INTERNAL;
+ }
+ else
+ {
+ edgeStat[edgeI] = EXTERNAL;
+ }
+ }
+ }
+ }
+ }
+}
+
+
// Returns next feature edge connected to pointI with correct value.
Foam::label Foam::surfaceFeatures::nextFeatEdge
(
@@ -436,6 +497,80 @@ Foam::surfaceFeatures::surfaceFeatures
}
+Foam::surfaceFeatures::surfaceFeatures
+(
+ const triSurface& surf,
+ const pointField& points,
+ const edgeList& edges,
+ const scalar mergeTol
+)
+:
+ surf_(surf),
+ featurePoints_(0),
+ featureEdges_(0),
+ externalStart_(0),
+ internalStart_(0)
+{
+ // Match edge mesh edges with the triSurface edges
+
+ const labelListList& surfEdgeFaces = surf_.edgeFaces();
+ const edgeList& surfEdges = surf_.edges();
+
+ scalar mergeTolSqr = sqr(mergeTol);
+
+ EdgeMap<label> dynFeatEdges(2*edges.size());
+ DynamicList<labelList> dynFeatureEdgeFaces(edges.size());
+
+ labelList edgeLabel;
+
+ nearestFeatEdge
+ (
+ edges,
+ points,
+ mergeTolSqr,
+ edgeLabel // label of surface edge or -1
+ );
+
+ label count = 0;
+ forAll(edgeLabel, sEdgeI)
+ {
+ const label sEdge = edgeLabel[sEdgeI];
+
+ if (sEdge == -1)
+ {
+ continue;
+ }
+
+ dynFeatEdges.insert(surfEdges[sEdge], count++);
+ dynFeatureEdgeFaces.append(surfEdgeFaces[sEdge]);
+ }
+
+ // Find whether an edge is external or internal
+ List<edgeStatus> edgeStat(dynFeatEdges.size(), NONE);
+
+ classifyFeatureAngles(dynFeatureEdgeFaces, edgeStat, GREAT);
+
+ // Transfer the edge status to a list encompassing all edges in the surface
+ // so that calcFeatPoints can be used.
+ List<edgeStatus> allEdgeStat(surf_.nEdges(), NONE);
+
+ forAll(allEdgeStat, eI)
+ {
+ EdgeMap<label>::const_iterator iter = dynFeatEdges.find(surfEdges[eI]);
+
+ if (iter != dynFeatEdges.end())
+ {
+ allEdgeStat[eI] = edgeStat[iter()];
+ }
+ }
+
+ edgeStat.clear();
+ dynFeatEdges.clear();
+
+ setFromStatus(allEdgeStat);
+}
+
+
//- Construct as copy
Foam::surfaceFeatures::surfaceFeatures(const surfaceFeatures& sf)
:
@@ -496,54 +631,10 @@ void Foam::surfaceFeatures::findFeatures(const scalar includedAngle)
{
scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
- const labelListList& edgeFaces = surf_.edgeFaces();
- const vectorField& faceNormals = surf_.faceNormals();
- const pointField& points = surf_.points();
-
// Per edge whether is feature edge.
List<edgeStatus> edgeStat(surf_.nEdges(), NONE);
- forAll(edgeFaces, edgeI)
- {
- const labelList& eFaces = edgeFaces[edgeI];
-
- if (eFaces.size() != 2)
- {
- // Non-manifold. What to do here? Is region edge? external edge?
- edgeStat[edgeI] = REGION;
- }
- else
- {
- label face0 = eFaces[0];
- label face1 = eFaces[1];
-
- if (surf_[face0].region() != surf_[face1].region())
- {
- edgeStat[edgeI] = REGION;
- }
- else
- {
- if ((faceNormals[face0] & faceNormals[face1]) < minCos)
- {
-
- // Check if convex or concave by looking at angle
- // between face centres and normal
- vector f0Tof1 =
- surf_[face1].centre(points)
- - surf_[face0].centre(points);
-
- if ((f0Tof1 & faceNormals[face0]) > 0.0)
- {
- edgeStat[edgeI] = INTERNAL;
- }
- else
- {
- edgeStat[edgeI] = EXTERNAL;
- }
- }
- }
- }
- }
+ classifyFeatureAngles(surf_.edgeFaces(), edgeStat, minCos);
setFromStatus(edgeStat);
}
@@ -1147,6 +1238,9 @@ void Foam::surfaceFeatures::nearestSurfEdge
const pointField& localPoints = surf_.localPoints();
+ treeBoundBox searchDomain(localPoints);
+ searchDomain.inflate(0.1);
+
indexedOctree<treeDataEdge> ppTree
(
treeDataEdge
@@ -1156,7 +1250,7 @@ void Foam::surfaceFeatures::nearestSurfEdge
localPoints,
selectedEdges
), // all information needed to do geometric checks
- treeBoundBox(localPoints), // overall search domain
+ searchDomain, // overall search domain
8, // maxLevel
10, // leafsize
3.0 // duplicity
@@ -1216,6 +1310,8 @@ void Foam::surfaceFeatures::nearestSurfEdge
pointOnEdge.setSize(selectedSampleEdges.size());
pointOnFeature.setSize(selectedSampleEdges.size());
+ treeBoundBox searchDomain(surf_.localPoints());
+
indexedOctree<treeDataEdge> ppTree
(
treeDataEdge
@@ -1224,11 +1320,11 @@ void Foam::surfaceFeatures::nearestSurfEdge
surf_.edges(),
surf_.localPoints(),
selectedEdges
- ), // all information needed to do geometric checks
- treeBoundBox(surf_.localPoints()), // overall search domain
- 8, // maxLevel
- 10, // leafsize
- 3.0 // duplicity
+ ), // all information needed to do geometric checks
+ searchDomain, // overall search domain
+ 8, // maxLevel
+ 10, // leafsize
+ 3.0 // duplicity
);
forAll(selectedSampleEdges, i)
@@ -1262,6 +1358,67 @@ void Foam::surfaceFeatures::nearestSurfEdge
}
+void Foam::surfaceFeatures::nearestFeatEdge
+(
+ const edgeList& edges,
+ const pointField& points,
+ scalar searchSpanSqr, // Search span
+ labelList& edgeLabel
+) const
+{
+ edgeLabel = labelList(surf_.nEdges(), -1);
+
+ treeBoundBox searchDomain(points);
+ searchDomain.inflate(0.1);
+
+ indexedOctree<treeDataEdge> ppTree
+ (
+ treeDataEdge
+ (
+ false,
+ edges,
+ points,
+ identity(edges.size())
+ ), // all information needed to do geometric checks
+ searchDomain, // overall search domain
+ 8, // maxLevel
+ 10, // leafsize
+ 3.0 // duplicity
+ );
+
+ const edgeList& surfEdges = surf_.edges();
+ const pointField& surfLocalPoints = surf_.localPoints();
+
+ forAll(surfEdges, edgeI)
+ {
+ const edge& sample = surfEdges[edgeI];
+
+ const point& startPoint = surfLocalPoints[sample.start()];
+ const point& midPoint = sample.centre(surfLocalPoints);
+
+ pointIndexHit infoMid = ppTree.findNearest
+ (
+ midPoint,
+ searchSpanSqr
+ );
+
+ if (infoMid.hit())
+ {
+ const vector surfEdgeDir = midPoint - startPoint;
+
+ const edge& featEdge = edges[infoMid.index()];
+ const vector featEdgeDir = featEdge.vec(points);
+
+ // Check that the edges are nearly parallel
+ if (mag(surfEdgeDir ^ featEdgeDir) < parallelTolerance)
+ {
+ edgeLabel[edgeI] = edgeI;
+ }
+ }
+ }
+}
+
+
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void Foam::surfaceFeatures::operator=(const surfaceFeatures& rhs)
diff --git a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H
index 7f7101b662a04046d8888d0b698f09735812141d..0fca9798a2e312c613ace576f1ba12b19c443069 100644
--- a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H
+++ b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H
@@ -95,6 +95,11 @@ private:
{}
};
+ // Static data
+
+ //- Tolerance for determining whether two vectors are parallel
+ static const scalar parallelTolerance;
+
// Private data
@@ -130,6 +135,14 @@ private:
//- Construct feature points where more than 2 feature edges meet
void calcFeatPoints(const List<edgeStatus>&);
+ //- Classify the angles of the feature edges
+ void classifyFeatureAngles
+ (
+ const labelListList& edgeFaces,
+ List<edgeStatus>& edgeStat,
+ const scalar minCos
+ ) const;
+
//- Choose next unset feature edge.
label nextFeatEdge
(
@@ -186,6 +199,15 @@ public:
//- Construct from file
surfaceFeatures(const triSurface&, const fileName& fName);
+ //- Construct from pointField and edgeList (edgeMesh)
+ surfaceFeatures
+ (
+ const triSurface&,
+ const pointField& points,
+ const edgeList& edges,
+ const scalar mergeTol = 1e-6
+ );
+
//- Construct as copy
surfaceFeatures(const surfaceFeatures&);
@@ -327,9 +349,8 @@ public:
pointField& edgePoint
) const;
-
//- Find nearest surface edge (out of selectedEdges) for each
- // sample edge.
+ // sample edge.
// Sets:
// - edgeLabel : label of surface edge.
// - pointOnEdge : exact position of nearest point on edge.
@@ -347,6 +368,16 @@ public:
pointField& pointOnFeature // point on sample edge
) const;
+ //- Find nearest feature edge to each surface edge. Uses the
+ // mid-point of the surface edges.
+ void nearestFeatEdge
+ (
+ const edgeList& edges,
+ const pointField& points,
+ scalar searchSpanSqr,
+ labelList& edgeLabel
+ ) const;
+
// Write