From 84fa233c009abc4d369c9fc6f7a1fe633257f8b1 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Mon, 22 Nov 2010 10:48:07 +0000
Subject: [PATCH] ENH: featureEdgeMesh : merged dev_cvm functionality.
---
.../surfaceFeatureExtract/Make/options | 3 +-
.../surfaceFeatureExtract.C | 23 +-
src/edgeMesh/Make/options | 6 +-
src/edgeMesh/edgeMesh.C | 10 -
src/edgeMesh/edgeMeshI.H | 10 +
.../featureEdgeMesh/featureEdgeMesh.C | 935 +++++++++++++++++-
.../featureEdgeMesh/featureEdgeMesh.H | 292 +++++-
7 files changed, 1234 insertions(+), 45 deletions(-)
diff --git a/applications/utilities/surface/surfaceFeatureExtract/Make/options b/applications/utilities/surface/surfaceFeatureExtract/Make/options
index 02e49d32c6e..f87dafaa5d1 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/Make/options
+++ b/applications/utilities/surface/surfaceFeatureExtract/Make/options
@@ -1,8 +1,9 @@
EXE_INC = \
- -I$(LIB_SRC)/cfdTools/general/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
+ -I$(LIB_SRC)/edgeMesh/lnInclude \
-I$(LIB_SRC)/triSurface/lnInclude
EXE_LIBS = \
-lmeshTools \
+ -ledgeMesh \
-ltriSurface
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index 83d84b92c82..7ba05f9d0e5 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -29,10 +29,13 @@ Description
\*---------------------------------------------------------------------------*/
+
#include "triangle.H"
#include "triSurface.H"
#include "argList.H"
+#include "Time.H"
#include "surfaceFeatures.H"
+#include "featureEdgeMesh.H"
#include "treeBoundBox.H"
#include "meshTools.H"
#include "OFstream.H"
@@ -135,7 +138,8 @@ int main(int argc, char *argv[])
"remove edges within specified bounding box"
);
- argList args(argc, argv);
+# include "setRootCase.H"
+# include "createTime.H"
Info<< "Feature line extraction is only valid on closed manifold surfaces."
<< endl;
@@ -282,7 +286,6 @@ int main(int argc, char *argv[])
Info<< endl << "Writing edge objs." << endl;
newSet.writeObj("final");
-
Info<< nl
<< "Final feature set:" << nl
<< " feature points : " << newSet.featurePoints().size() << nl
@@ -293,6 +296,22 @@ int main(int argc, char *argv[])
<< " internal edges : " << newSet.nInternalEdges() << nl
<< endl;
+ // Extracting and writing a featureEdgeMesh
+
+ Pout<< nl << "Writing featureEdgeMesh to constant/featureEdgeMesh."
+ << endl;
+
+ featureEdgeMesh feMesh
+ (
+ newSet,
+ runTime,
+ surfFileName.lessExt().name() + ".featureEdgeMesh"
+ );
+
+ feMesh.writeObj(surfFileName.lessExt().name());
+
+ feMesh.write();
+
Info<< "End\n" << endl;
return 0;
diff --git a/src/edgeMesh/Make/options b/src/edgeMesh/Make/options
index 7e207d0dbea..4796c18e22c 100644
--- a/src/edgeMesh/Make/options
+++ b/src/edgeMesh/Make/options
@@ -1,5 +1,9 @@
EXE_INC = \
- -I$(LIB_SRC)/fileFormats/lnInclude
+ -I$(LIB_SRC)/fileFormats/lnInclude \
+ -I$(LIB_SRC)/triSurface/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude
LIB_LIBS = \
+ -ltriSurface \
+ -lmeshTools \
-lfileFormats
diff --git a/src/edgeMesh/edgeMesh.C b/src/edgeMesh/edgeMesh.C
index d5c8f12ef8f..7dd70151445 100644
--- a/src/edgeMesh/edgeMesh.C
+++ b/src/edgeMesh/edgeMesh.C
@@ -387,14 +387,4 @@ void Foam::edgeMesh::mergePoints(const scalar mergeDist)
}
-// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
-
-void Foam::edgeMesh::operator=(const edgeMesh& rhs)
-{
- points_ = rhs.points_;
- edges_ = rhs.edges_;
- pointEdgesPtr_.clear();
-}
-
-
// ************************************************************************* //
diff --git a/src/edgeMesh/edgeMeshI.H b/src/edgeMesh/edgeMeshI.H
index d5e38343ec4..c15fbd14ce9 100644
--- a/src/edgeMesh/edgeMeshI.H
+++ b/src/edgeMesh/edgeMeshI.H
@@ -69,4 +69,14 @@ inline Foam::edgeList& Foam::edgeMesh::storedEdges()
}
+// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
+
+void Foam::edgeMesh::operator=(const edgeMesh& rhs)
+{
+ points_ = rhs.points_;
+ edges_ = rhs.edges_;
+ pointEdgesPtr_.clear();
+}
+
+
// ************************************************************************* //
diff --git a/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C
index 34bcd4388e4..da89dcf2d53 100644
--- a/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C
+++ b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C
@@ -24,18 +24,56 @@ License
\*---------------------------------------------------------------------------*/
#include "featureEdgeMesh.H"
+#include "triSurface.H"
+#include "Random.H"
+#include "Time.H"
+#include "meshTools.H"
+#include "linePointRef.H"
+#include "ListListOps.H"
+#include "OFstream.H"
+#include "IFstream.H"
+#include "unitConversion.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::featureEdgeMesh, 0);
+Foam::scalar Foam::featureEdgeMesh::cosNormalAngleTol_ =
+ Foam::cos(degToRad(0.1));
+
+
+Foam::label Foam::featureEdgeMesh::convexStart_ = 0;
+
+
+Foam::label Foam::featureEdgeMesh::externalStart_ = 0;
+
+
+Foam::label Foam::featureEdgeMesh::nPointTypes = 4;
+
+
+Foam::label Foam::featureEdgeMesh::nEdgeTypes = 5;
+
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::featureEdgeMesh::featureEdgeMesh(const IOobject& io)
:
regIOobject(io),
- edgeMesh(pointField(0), edgeList(0))
+ edgeMesh(pointField(0), edgeList(0)),
+ concaveStart_(0),
+ mixedStart_(0),
+ nonFeatureStart_(0),
+ internalStart_(0),
+ flatStart_(0),
+ openStart_(0),
+ multipleStart_(0),
+ normals_(0),
+ edgeDirections_(0),
+ edgeNormals_(0),
+ featurePointNormals_(0),
+ regionEdges_(0),
+ edgeTree_(),
+ edgeTreesByType_()
{
if
(
@@ -52,9 +90,39 @@ Foam::featureEdgeMesh::featureEdgeMesh(const IOobject& io)
<< endl;
}
+ Istream& is = readStream(typeName);
+
+ is >> *this
+ >> concaveStart_
+ >> mixedStart_
+ >> nonFeatureStart_
+ >> internalStart_
+ >> flatStart_
+ >> openStart_
+ >> multipleStart_
+ >> normals_
+ >> edgeNormals_
+ >> featurePointNormals_
+ >> regionEdges_;
- readStream(typeName) >> *this;
close();
+
+ {
+ // Calculate edgeDirections
+
+ const edgeList& eds(edges());
+
+ const pointField& pts(points());
+
+ edgeDirections_.setSize(eds.size());
+
+ forAll(eds, eI)
+ {
+ edgeDirections_[eI] = eds[eI].vec(pts);
+ }
+
+ edgeDirections_ /= mag(edgeDirections_);
+ }
}
if (debug)
@@ -71,11 +139,25 @@ Foam::featureEdgeMesh::featureEdgeMesh(const IOobject& io)
Foam::featureEdgeMesh::featureEdgeMesh
(
const IOobject& io,
- const featureEdgeMesh& em
+ const featureEdgeMesh& fem
)
:
regIOobject(io),
- edgeMesh(em)
+ edgeMesh(fem),
+ concaveStart_(fem.concaveStart()),
+ mixedStart_(fem.mixedStart()),
+ nonFeatureStart_(fem.nonFeatureStart()),
+ internalStart_(fem.internalStart()),
+ flatStart_(fem.flatStart()),
+ openStart_(fem.openStart()),
+ multipleStart_(fem.multipleStart()),
+ normals_(fem.normals()),
+ edgeDirections_(fem.edgeDirections()),
+ edgeNormals_(fem.edgeNormals()),
+ featurePointNormals_(fem.featurePointNormals()),
+ regionEdges_(fem.regionEdges()),
+ edgeTree_(),
+ edgeTreesByType_()
{}
@@ -87,42 +169,851 @@ Foam::featureEdgeMesh::featureEdgeMesh
)
:
regIOobject(io),
- edgeMesh(pointLst, edgeLst)
+ edgeMesh(pointLst, edgeLst),
+ concaveStart_(0),
+ mixedStart_(0),
+ nonFeatureStart_(0),
+ internalStart_(0),
+ flatStart_(0),
+ openStart_(0),
+ multipleStart_(0),
+ normals_(0),
+ edgeDirections_(0),
+ edgeNormals_(0),
+ featurePointNormals_(0),
+ regionEdges_(0),
+ edgeTree_(),
+ edgeTreesByType_()
+{}
+
+
+Foam::featureEdgeMesh::featureEdgeMesh
+(
+ const surfaceFeatures& sFeat,
+ const objectRegistry& obr,
+ const fileName& sFeatFileName
+)
+:
+ regIOobject
+ (
+ IOobject
+ (
+ sFeatFileName,
+ obr.time().constant(),
+ "featureEdgeMesh",
+ obr,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ )
+ ),
+ edgeMesh(pointField(0), edgeList(0)),
+ concaveStart_(-1),
+ mixedStart_(-1),
+ nonFeatureStart_(-1),
+ internalStart_(-1),
+ flatStart_(-1),
+ openStart_(-1),
+ multipleStart_(-1),
+ normals_(0),
+ edgeDirections_(0),
+ edgeNormals_(0),
+ featurePointNormals_(0),
+ regionEdges_(0),
+ edgeTree_(),
+ edgeTreesByType_()
{
- if
+ // Extract and reorder the data from surfaceFeatures
+
+ // References to the surfaceFeatures data
+ const triSurface& surf(sFeat.surface());
+ const pointField& sFeatLocalPts(surf.localPoints());
+ const edgeList& sFeatEds(surf.edges());
+
+ // Filling the featureEdgeMesh with the raw geometrical data.
+
+ label nFeatEds = sFeat.featureEdges().size();
+
+ DynamicList<point> tmpPts;
+ edgeList eds(nFeatEds);
+ DynamicList<vector> norms;
+ vectorField edgeDirections(nFeatEds);
+ labelListList edgeNormals(nFeatEds);
+ DynamicList<label> regionEdges;
+
+ // Mapping between old and new indices, there is entry in the map for each
+ // of sFeatLocalPts, -1 means that this point hasn't been used (yet), >= 0
+ // corresponds to the index
+ labelList pointMap(sFeatLocalPts.size(), -1);
+
+ // Noting when the normal of a face has been used so not to duplicate
+ labelList faceMap(surf.size(), -1);
+
+ // Collecting the status of edge for subsequent sorting
+ List<edgeStatus> edStatus(nFeatEds, NONE);
+
+ forAll(sFeat.featurePoints(), i)
+ {
+ label sFPI = sFeat.featurePoints()[i];
+
+ tmpPts.append(sFeatLocalPts[sFPI]);
+
+ pointMap[sFPI] = tmpPts.size() - 1;
+ }
+
+ // All feature points have been added
+ nonFeatureStart_ = tmpPts.size();
+
+ forAll(sFeat.featureEdges(), i)
+ {
+ label sFEI = sFeat.featureEdges()[i];
+
+ const edge& fE(sFeatEds[sFEI]);
+
+ // Check to see if the points have been already used
+ if (pointMap[fE.start()] == -1)
+ {
+ tmpPts.append(sFeatLocalPts[fE.start()]);
+
+ pointMap[fE.start()] = tmpPts.size() - 1;
+ }
+
+ eds[i].start() = pointMap[fE.start()];
+
+ if (pointMap[fE.end()] == -1)
+ {
+ tmpPts.append(sFeatLocalPts[fE.end()]);
+
+ pointMap[fE.end()] = tmpPts.size() - 1;
+ }
+
+ eds[i].end() = pointMap[fE.end()];
+
+ // Pick up the faces adjacent to the feature edge
+ const labelList& eFaces = surf.edgeFaces()[sFEI];
+
+ edgeNormals[i].setSize(eFaces.size());
+
+ forAll(eFaces, j)
+ {
+ label eFI = eFaces[j];
+
+ // Check to see if the points have been already used
+ if (faceMap[eFI] == -1)
+ {
+ norms.append(surf.faceNormals()[eFI]);
+
+ faceMap[eFI] = norms.size() - 1;
+ }
+
+ edgeNormals[i][j] = faceMap[eFI];
+ }
+
+ vector fC0tofC1(vector::zero);
+
+ if (eFaces.size() == 2)
+ {
+ fC0tofC1 =
+ surf[eFaces[1]].centre(surf.points())
+ - surf[eFaces[0]].centre(surf.points());
+ }
+
+ edStatus[i] = classifyEdge(norms, edgeNormals[i], fC0tofC1);
+
+ edgeDirections[i] = fE.vec(sFeatLocalPts);
+
+ if (i < sFeat.nRegionEdges())
+ {
+ regionEdges.append(i);
+ }
+ }
+
+ // Reorder the edges by classification
+
+ List<DynamicList<label> > allEds(nEdgeTypes);
+
+ DynamicList<label>& externalEds(allEds[0]);
+ DynamicList<label>& internalEds(allEds[1]);
+ DynamicList<label>& flatEds(allEds[2]);
+ DynamicList<label>& openEds(allEds[3]);
+ DynamicList<label>& multipleEds(allEds[4]);
+
+ forAll(eds, i)
+ {
+ edgeStatus eStat = edStatus[i];
+
+ if (eStat == EXTERNAL)
+ {
+ externalEds.append(i);
+ }
+ else if (eStat == INTERNAL)
+ {
+ internalEds.append(i);
+ }
+ else if (eStat == FLAT)
+ {
+ flatEds.append(i);
+ }
+ else if (eStat == OPEN)
+ {
+ openEds.append(i);
+ }
+ else if (eStat == MULTIPLE)
+ {
+ multipleEds.append(i);
+ }
+ else if (eStat == NONE)
+ {
+ FatalErrorIn
+ (
+ "Foam::featureEdgeMesh::featureEdgeMesh"
+ "("
+ " const surfaceFeatures& sFeat,"
+ " const objectRegistry& obr,"
+ " const fileName& sFeatFileName"
+ ")"
+ )
+ << nl << "classifyEdge returned NONE on edge "
+ << eds[i]
+ << ". There is a problem with definition of this edge."
+ << nl << abort(FatalError);
+ }
+ }
+
+ internalStart_ = externalEds.size();
+ flatStart_ = internalStart_ + internalEds.size();
+ openStart_ = flatStart_ + flatEds.size();
+ multipleStart_ = openStart_ + openEds.size();
+
+ labelList edMap
(
- io.readOpt() == IOobject::MUST_READ
- || io.readOpt() == IOobject::MUST_READ_IF_MODIFIED
- || (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
- )
+ ListListOps::combine<labelList>
+ (
+ allEds,
+ accessOp<labelList>()
+ )
+ );
+
+ edMap = invert(edMap.size(), edMap);
+
+ inplaceReorder(edMap, eds);
+ inplaceReorder(edMap, edStatus);
+ inplaceReorder(edMap, edgeDirections);
+ inplaceReorder(edMap, edgeNormals);
+ inplaceRenumber(edMap, regionEdges);
+
+ pointField pts(tmpPts);
+
+ // Initialise the edgeMesh
+ edgeMesh::operator=(edgeMesh(pts, eds));
+
+ // Initialise sorted edge related data
+ edgeDirections_ = edgeDirections/mag(edgeDirections);
+ edgeNormals_ = edgeNormals;
+ regionEdges_ = regionEdges;
+
+ // Normals are all now found and indirectly addressed, can also be stored
+ normals_ = vectorField(norms);
+
+ // Reorder the feature points by classification
+
+ List<DynamicList<label> > allPts(3);
+
+ DynamicList<label>& convexPts(allPts[0]);
+ DynamicList<label>& concavePts(allPts[1]);
+ DynamicList<label>& mixedPts(allPts[2]);
+
+ for (label i = 0; i < nonFeatureStart_; i++)
{
- if (readOpt() == IOobject::MUST_READ_IF_MODIFIED)
+ pointStatus ptStatus = classifyFeaturePoint(i);
+
+ if (ptStatus == CONVEX)
{
- WarningIn("featureEdgeMesh::featureEdgeMesh(const IOobject&)")
- << "Specified IOobject::MUST_READ_IF_MODIFIED but class"
- << " does not support automatic rereading."
- << endl;
+ convexPts.append(i);
+ }
+ else if (ptStatus == CONCAVE)
+ {
+ concavePts.append(i);
}
+ else if (ptStatus == MIXED)
+ {
+ mixedPts.append(i);
+ }
+ else if (ptStatus == NONFEATURE)
+ {
+ FatalErrorIn
+ (
+ "Foam::featureEdgeMesh::featureEdgeMesh"
+ "("
+ " const surfaceFeatures& sFeat,"
+ " const objectRegistry& obr,"
+ " const fileName& sFeatFileName"
+ ")"
+ )
+ << nl << "classifyFeaturePoint returned NONFEATURE on point at "
+ << points()[i]
+ << ". There is a problem with definition of this feature point."
+ << nl << abort(FatalError);
+ }
+ }
+ concaveStart_ = convexPts.size();
+ mixedStart_ = concaveStart_ + concavePts.size();
- readStream(typeName) >> *this;
- close();
+ labelList ftPtMap
+ (
+ ListListOps::combine<labelList>
+ (
+ allPts,
+ accessOp<labelList>()
+ )
+ );
+
+ ftPtMap = invert(ftPtMap.size(), ftPtMap);
+
+ // Creating the ptMap from the ftPtMap with identity values up to the size
+ // of pts to create an oldToNew map for inplaceReorder
+
+ labelList ptMap(identity(pts.size()));
+
+ forAll(ftPtMap, i)
+ {
+ ptMap[i] = ftPtMap[i];
+ }
+
+ inplaceReorder(ptMap, pts);
+
+ forAll(eds, i)
+ {
+ inplaceRenumber(ptMap, eds[i]);
+ }
+
+ // Reinitialise the edgeMesh with sorted feature points and
+ // renumbered edges
+ edgeMesh::operator=(edgeMesh(pts, eds));
+
+ // Generate the featurePointNormals
+
+ labelListList featurePointNormals(nonFeatureStart_);
+
+ for (label i = 0; i < nonFeatureStart_; i++)
+ {
+ DynamicList<label> tmpFtPtNorms;
+
+ const labelList& ptEds = pointEdges()[i];
+
+ forAll(ptEds, j)
+ {
+ const labelList& ptEdNorms(edgeNormals[ptEds[j]]);
+
+ forAll(ptEdNorms, k)
+ {
+ if (findIndex(tmpFtPtNorms, ptEdNorms[k]) == -1)
+ {
+ bool addNormal = true;
+
+ // Check that the normal direction is unique at this feature
+ forAll(tmpFtPtNorms, q)
+ {
+ if
+ (
+ (normals_[ptEdNorms[k]] & normals_[tmpFtPtNorms[q]])
+ > cosNormalAngleTol_
+ )
+ {
+ // Parallel to an existing normal, do not add
+ addNormal = false;
+
+ break;
+ }
+ }
+
+ if (addNormal)
+ {
+ tmpFtPtNorms.append(ptEdNorms[k]);
+ }
+ }
+ }
+ }
+
+ featurePointNormals[i] = tmpFtPtNorms;
+ }
+
+ featurePointNormals_ = featurePointNormals;
+}
+
+
+Foam::featureEdgeMesh::featureEdgeMesh
+(
+ const IOobject& io,
+ const pointField& pts,
+ const edgeList& eds,
+ label concaveStart,
+ label mixedStart,
+ label nonFeatureStart,
+ label internalStart,
+ label flatStart,
+ label openStart,
+ label multipleStart,
+ const vectorField& normals,
+ const vectorField& edgeDirections,
+ const labelListList& edgeNormals,
+ const labelListList& featurePointNormals,
+ const labelList& regionEdges
+)
+:
+ regIOobject(io),
+ edgeMesh(pts, eds),
+ concaveStart_(concaveStart),
+ mixedStart_(mixedStart),
+ nonFeatureStart_(nonFeatureStart),
+ internalStart_(internalStart),
+ flatStart_(flatStart),
+ openStart_(openStart),
+ multipleStart_(multipleStart),
+ normals_(normals),
+ edgeDirections_(edgeDirections),
+ edgeNormals_(edgeNormals),
+ featurePointNormals_(featurePointNormals),
+ regionEdges_(regionEdges),
+ edgeTree_(),
+ edgeTreesByType_()
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::featureEdgeMesh::~featureEdgeMesh()
+{}
+
+
+// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
+
+Foam::featureEdgeMesh::pointStatus Foam::featureEdgeMesh::classifyFeaturePoint
+(
+ label ptI
+) const
+{
+ labelList ptEds(pointEdges()[ptI]);
+
+ label nPtEds = ptEds.size();
+ label nExternal = 0;
+ label nInternal = 0;
+
+ if (nPtEds == 0)
+ {
+ // There are no edges attached to the point, this is a problem
+ return NONFEATURE;
+ }
+
+ forAll(ptEds, i)
+ {
+ edgeStatus edStat = getEdgeStatus(ptEds[i]);
+
+ if (edStat == EXTERNAL)
+ {
+ nExternal++;
+ }
+ else if (edStat == INTERNAL)
+ {
+ nInternal++;
+ }
+ }
+
+ if (nExternal == nPtEds)
+ {
+ return CONVEX;
+ }
+ else if (nInternal == nPtEds)
+ {
+ return CONCAVE;
+ }
+ else
+ {
+ return MIXED;
+ }
+}
+
+
+Foam::featureEdgeMesh::edgeStatus Foam::featureEdgeMesh::classifyEdge
+(
+ const List<vector>& norms,
+ const labelList& edNorms,
+ const vector& fC0tofC1
+) const
+{
+ label nEdNorms = edNorms.size();
+
+ if (nEdNorms == 1)
+ {
+ return OPEN;
+ }
+ else if (nEdNorms == 2)
+ {
+ const vector n0(norms[edNorms[0]]);
+ const vector n1(norms[edNorms[1]]);
+
+ if ((n0 & n1) > cosNormalAngleTol_)
+ {
+ return FLAT;
+ }
+ else if ((fC0tofC1 & n0) > 0.0)
+ {
+ return INTERNAL;
+ }
+ else
+ {
+ return EXTERNAL;
+ }
+ }
+ else if (nEdNorms > 2)
+ {
+ return MULTIPLE;
+ }
+ else
+ {
+ // There is a problem - the edge has no normals
+ return NONE;
+ }
+}
+
+
+// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
+
+void Foam::featureEdgeMesh::nearestFeatureEdge
+(
+ const point& sample,
+ scalar searchDistSqr,
+ pointIndexHit& info
+) const
+{
+ info = edgeTree().findNearest
+ (
+ sample,
+ searchDistSqr
+ );
+}
+
+
+void Foam::featureEdgeMesh::nearestFeatureEdge
+(
+ const pointField& samples,
+ const scalarField& searchDistSqr,
+ List<pointIndexHit>& info
+) const
+{
+ info.setSize(samples.size());
+
+ forAll(samples, i)
+ {
+ nearestFeatureEdge
+ (
+ samples[i],
+ searchDistSqr[i],
+ info[i]
+ );
+ }
+}
+
+
+void Foam::featureEdgeMesh::nearestFeatureEdgeByType
+(
+ const point& sample,
+ const scalarField& searchDistSqr,
+ List<pointIndexHit>& info
+) const
+{
+ const PtrList<indexedOctree<treeDataEdge> >& edgeTrees = edgeTreesByType();
+
+ info.setSize(edgeTrees.size());
+
+ labelList sliceStarts(edgeTrees.size());
+
+ sliceStarts[0] = externalStart_;
+ sliceStarts[1] = internalStart_;
+ sliceStarts[2] = flatStart_;
+ sliceStarts[3] = openStart_;
+ sliceStarts[4] = multipleStart_;
+
+ forAll(edgeTrees, i)
+ {
+ info[i] = edgeTrees[i].findNearest
+ (
+ sample,
+ searchDistSqr[i]
+ );
+
+ // The index returned by the indexedOctree is local to the slice of
+ // edges it was supplied with, return the index to the value in the
+ // complete edge list
+
+ info[i].setIndex(info[i].index() + sliceStarts[i]);
+ }
+}
+
+
+const Foam::indexedOctree<Foam::treeDataEdge>&
+Foam::featureEdgeMesh::edgeTree() const
+{
+ if (edgeTree_.empty())
+ {
+ Random rndGen(17301893);
+
+ // Slightly extended bb. Slightly off-centred just so on symmetric
+ // geometry there are less face/edge aligned items.
+ treeBoundBox bb
+ (
+ treeBoundBox(points()).extend(rndGen, 1E-4)
+ );
+
+ bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+ bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+ labelList allEdges(identity(edges().size()));
+
+ edgeTree_.reset
+ (
+ new indexedOctree<treeDataEdge>
+ (
+ treeDataEdge
+ (
+ false, // cachebb
+ edges(), // edges
+ points(), // points
+ allEdges // selected edges
+ ),
+ bb, // bb
+ 8, // maxLevel
+ 10, // leafsize
+ 3.0 // duplicity
+ )
+ );
}
+
+ return edgeTree_();
}
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+const Foam::PtrList<Foam::indexedOctree<Foam::treeDataEdge> >&
+Foam::featureEdgeMesh::edgeTreesByType() const
+{
+ if (edgeTreesByType_.size() == 0)
+ {
+ edgeTreesByType_.setSize(nEdgeTypes);
+
+ Random rndGen(872141);
+
+ // Slightly extended bb. Slightly off-centred just so on symmetric
+ // geometry there are less face/edge aligned items.
+ treeBoundBox bb
+ (
+ treeBoundBox(points()).extend(rndGen, 1E-4)
+ );
+
+ bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+ bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+ labelListList sliceEdges(nEdgeTypes);
+
+ // External edges
+ sliceEdges[0] =
+ identity(internalStart_ - externalStart_) + externalStart_;
+
+ // Internal edges
+ sliceEdges[1] = identity(flatStart_ - internalStart_) + internalStart_;
+
+ // Flat edges
+ sliceEdges[2] = identity(openStart_ - flatStart_) + flatStart_;
+
+ // Open edges
+ sliceEdges[3] = identity(multipleStart_ - openStart_) + openStart_;
+
+ // Multiple edges
+ sliceEdges[4] =
+ identity(edges().size() - multipleStart_) + multipleStart_;
+
+ forAll(edgeTreesByType_, i)
+ {
+ edgeTreesByType_.set
+ (
+ i,
+ new indexedOctree<treeDataEdge>
+ (
+ treeDataEdge
+ (
+ false, // cachebb
+ edges(), // edges
+ points(), // points
+ sliceEdges[i] // selected edges
+ ),
+ bb, // bb
+ 8, // maxLevel
+ 10, // leafsize
+ 3.0 // duplicity
+ )
+ );
+ }
+ }
+
+ return edgeTreesByType_;
+}
+
-bool Foam::featureEdgeMesh::readData(Istream& is)
+void Foam::featureEdgeMesh::writeObj
+(
+ const fileName& prefix
+) const
{
- is >> *this;
- return !is.bad();
+ Pout<< nl << "Writing featureEdgeMesh components to " << prefix << endl;
+
+ label verti = 0;
+
+ edgeMesh::write(prefix + "_edgeMesh.obj");
+
+ OFstream convexFtPtStr(prefix + "_convexFeaturePts.obj");
+ Pout<< "Writing convex feature points to " << convexFtPtStr.name() << endl;
+
+ for(label i = 0; i < concaveStart_; i++)
+ {
+ meshTools::writeOBJ(convexFtPtStr, points()[i]);
+ }
+
+ OFstream concaveFtPtStr(prefix + "_concaveFeaturePts.obj");
+ Pout<< "Writing concave feature points to "
+ << concaveFtPtStr.name() << endl;
+
+ for(label i = concaveStart_; i < mixedStart_; i++)
+ {
+ meshTools::writeOBJ(concaveFtPtStr, points()[i]);
+ }
+
+ OFstream mixedFtPtStr(prefix + "_mixedFeaturePts.obj");
+ Pout<< "Writing mixed feature points to " << mixedFtPtStr.name() << endl;
+
+ for(label i = mixedStart_; i < nonFeatureStart_; i++)
+ {
+ meshTools::writeOBJ(mixedFtPtStr, points()[i]);
+ }
+
+ OFstream mixedFtPtStructureStr(prefix + "_mixedFeaturePtsStructure.obj");
+ Pout<< "Writing mixed feature point structure to "
+ << mixedFtPtStructureStr.name() << endl;
+
+ verti = 0;
+ for(label i = mixedStart_; i < nonFeatureStart_; i++)
+ {
+ const labelList& ptEds = pointEdges()[i];
+
+ forAll(ptEds, j)
+ {
+ const edge& e = edges()[ptEds[j]];
+
+ meshTools::writeOBJ(mixedFtPtStructureStr, points()[e[0]]); verti++;
+ meshTools::writeOBJ(mixedFtPtStructureStr, points()[e[1]]); verti++;
+ mixedFtPtStructureStr << "l " << verti-1 << ' ' << verti << endl;
+ }
+ }
+
+ OFstream externalStr(prefix + "_externalEdges.obj");
+ Pout<< "Writing external edges to " << externalStr.name() << endl;
+
+ verti = 0;
+ for (label i = externalStart_; i < internalStart_; i++)
+ {
+ const edge& e = edges()[i];
+
+ meshTools::writeOBJ(externalStr, points()[e[0]]); verti++;
+ meshTools::writeOBJ(externalStr, points()[e[1]]); verti++;
+ externalStr << "l " << verti-1 << ' ' << verti << endl;
+ }
+
+ OFstream internalStr(prefix + "_internalEdges.obj");
+ Pout<< "Writing internal edges to " << internalStr.name() << endl;
+
+ verti = 0;
+ for (label i = internalStart_; i < flatStart_; i++)
+ {
+ const edge& e = edges()[i];
+
+ meshTools::writeOBJ(internalStr, points()[e[0]]); verti++;
+ meshTools::writeOBJ(internalStr, points()[e[1]]); verti++;
+ internalStr << "l " << verti-1 << ' ' << verti << endl;
+ }
+
+ OFstream flatStr(prefix + "_flatEdges.obj");
+ Pout<< "Writing flat edges to " << flatStr.name() << endl;
+
+ verti = 0;
+ for (label i = flatStart_; i < openStart_; i++)
+ {
+ const edge& e = edges()[i];
+
+ meshTools::writeOBJ(flatStr, points()[e[0]]); verti++;
+ meshTools::writeOBJ(flatStr, points()[e[1]]); verti++;
+ flatStr << "l " << verti-1 << ' ' << verti << endl;
+ }
+
+ OFstream openStr(prefix + "_openEdges.obj");
+ Pout<< "Writing open edges to " << openStr.name() << endl;
+
+ verti = 0;
+ for (label i = openStart_; i < multipleStart_; i++)
+ {
+ const edge& e = edges()[i];
+
+ meshTools::writeOBJ(openStr, points()[e[0]]); verti++;
+ meshTools::writeOBJ(openStr, points()[e[1]]); verti++;
+ openStr << "l " << verti-1 << ' ' << verti << endl;
+ }
+
+ OFstream multipleStr(prefix + "_multipleEdges.obj");
+ Pout<< "Writing multiple edges to " << multipleStr.name() << endl;
+
+ verti = 0;
+ for (label i = multipleStart_; i < edges().size(); i++)
+ {
+ const edge& e = edges()[i];
+
+ meshTools::writeOBJ(multipleStr, points()[e[0]]); verti++;
+ meshTools::writeOBJ(multipleStr, points()[e[1]]); verti++;
+ multipleStr << "l " << verti-1 << ' ' << verti << endl;
+ }
+
+ OFstream regionStr(prefix + "_regionEdges.obj");
+ Pout<< "Writing region edges to " << regionStr.name() << endl;
+
+ verti = 0;
+ forAll(regionEdges_, i)
+ {
+ const edge& e = edges()[regionEdges_[i]];
+
+ meshTools::writeOBJ(regionStr, points()[e[0]]); verti++;
+ meshTools::writeOBJ(regionStr, points()[e[1]]); verti++;
+ regionStr << "l " << verti-1 << ' ' << verti << endl;
+ }
}
bool Foam::featureEdgeMesh::writeData(Ostream& os) const
{
- os << *this;
+ os << "// points, edges, concaveStart, mixedStart, nonFeatureStart, " << nl
+ << "// internalStart, flatStart, openStart, multipleStart, " << nl
+ << "// normals, edgeNormals, featurePointNormals, regionEdges" << nl
+ << endl;
+
+ os << points() << nl
+ << edges() << nl
+ << concaveStart_ << token::SPACE
+ << mixedStart_ << token::SPACE
+ << nonFeatureStart_ << token::SPACE
+ << internalStart_ << token::SPACE
+ << flatStart_ << token::SPACE
+ << openStart_ << token::SPACE
+ << multipleStart_ << nl
+ << normals_ << nl
+ << edgeNormals_ << nl
+ << featurePointNormals_ << nl
+ << regionEdges_
+ << endl;
return os.good();
}
diff --git a/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H
index 613c9bdd94b..93f034d9aaf 100644
--- a/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H
+++ b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H
@@ -25,9 +25,27 @@ Class
Foam::featureEdgeMesh
Description
- features (lines), readable from file
+
+ Feature points are a sorted subset at the start of the overall points list:
+ 0 .. concaveStart_-1 : convex points (w.r.t normals)
+ concaveStart_-1 .. mixedStart_-1 : concave points
+ mixedStart_ .. nonFeatureStart_ : mixed internal/external points
+ nonFeatureStart_ .. size-1 : non-feature points
+
+ Feature edges are the edgeList of the edgeMesh and are sorted:
+ 0 .. internalStart_-1 : external edges (convex w.r.t normals)
+ internalStart_ .. flatStart_-1 : internal edges (concave)
+ flatStart_ .. openStart_-1 : flat edges (neither concave or convex)
+ can arise from region interfaces on
+ flat surfaces
+ openStart_ .. multipleStart_-1 : open edges (e.g. from baffle surfaces)
+ multipleStart_ .. size-1 : multiply connected edges
+
+ The edge direction and feature edge and feature point adjacent normals
+ are stored.
SourceFiles
+ featureEdgeMeshI.H
featureEdgeMesh.C
\*---------------------------------------------------------------------------*/
@@ -36,7 +54,12 @@ SourceFiles
#define featureEdgeMesh_H
#include "edgeMesh.H"
-#include "regIOobject.H"
+#include "surfaceFeatures.H"
+#include "objectRegistry.H"
+#include "IOdictionary.H"
+#include "indexedOctree.H"
+#include "treeDataEdge.H"
+#include "pointIndexHit.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -58,11 +81,115 @@ public:
//- Runtime type information
TypeName("featureEdgeMesh");
+ enum pointStatus
+ {
+ CONVEX, // Fully convex point (w.r.t normals)
+ CONCAVE, // Fully concave point
+ MIXED, // A point surrounded by both convex and concave edges
+ NONFEATURE // Not a feature point
+ };
+
+ enum edgeStatus
+ {
+ EXTERNAL, // "Convex" edge
+ INTERNAL, // "Concave" edge
+ FLAT, // Neither concave or convex, on a flat surface
+ OPEN, // i.e. only connected to one face
+ MULTIPLE, // Multiply connected (connected to more than two faces)
+ NONE // Not a classified feature edge (consistency with
+ // surfaceFeatures)
+ };
+
+private:
+
+ // Static data
+
+ //- Angular closeness tolerance for treating normals as the same
+ static scalar cosNormalAngleTol_;
+
+ //- Index of the start of the convex feature points - static as 0
+ static label convexStart_;
+
+ //- Index of the start of the external feature edges - static as 0
+ static label externalStart_;
+
+
+ // Private data
+
+ //- Index of the start of the concave feature points
+ label concaveStart_;
+
+ //- Index of the start of the mixed type feature points
+ label mixedStart_;
+
+ //- Index of the start of the non-feature points
+ label nonFeatureStart_;
+
+ //- Index of the start of the internal feature edges
+ label internalStart_;
+
+ //- Index of the start of the flat feature edges
+ label flatStart_;
+
+ //- Index of the start of the open feature edges
+ label openStart_;
+
+ //- Index of the start of the multiply-connected feature edges
+ label multipleStart_;
+
+ //- Normals of the features, to be referred to by index by both feature
+ // points and edges, unsorted
+ vectorField normals_;
+
+ //- Flat and open edges require the direction of the edge
+ vectorField edgeDirections_;
+
+ //- Indices of the normals that are adjacent to the feature edges
+ labelListList edgeNormals_;
+
+ //- Indices of the normals that are adjacent to the feature points
+ labelListList featurePointNormals_;
+
+ //- Feature edges which are on the boundary between regions
+ labelList regionEdges_;
+
+ //- Search tree for all edges
+ mutable autoPtr<indexedOctree<treeDataEdge> > edgeTree_;
+
+ //- Individual search trees for each type of edge
+ mutable PtrList<indexedOctree<treeDataEdge> > edgeTreesByType_;
+
+
+ // Private Member Functions
+
+ //- Classify the type of feature point. Requires valid stored member
+ // data for edges and normals.
+ pointStatus classifyFeaturePoint(label ptI) const;
+
+ //- Classify the type of feature edge. Requires face centre 0 to face
+ // centre 1 vector to distinguish internal from external
+ edgeStatus classifyEdge
+ (
+ const List<vector>& norms,
+ const labelList& edNorms,
+ const vector& fC0tofC1
+ ) const;
+
+
+public:
+
+ // Static data
+
+ //- Number of possible point types (i.e. number of slices)
+ static label nPointTypes;
+
+ //- Number of possible feature edge types (i.e. number of slices)
+ static label nEdgeTypes;
// Constructors
//- Construct (read) given an IOobject
- explicit featureEdgeMesh(const IOobject&);
+ featureEdgeMesh(const IOobject&);
//- Construct as copy
explicit featureEdgeMesh(const IOobject&, const featureEdgeMesh&);
@@ -75,14 +202,157 @@ public:
const Xfer<edgeList>&
);
- //- Give precedence to the regIOobject write
- using regIOobject::write;
+ //- Construct (read) given surfaceFeatures, an objectRegistry and a
+ // fileName to write to. Extracts, classifies and reorders the data
+ // from surfaceFeatures.
+ featureEdgeMesh
+ (
+ const surfaceFeatures& sFeat,
+ const objectRegistry& obr,
+ const fileName& sFeatFileName
+ );
- //- ReadData function required for regIOobject read operation
- virtual bool readData(Istream&);
+ //- Construct from all components
+ featureEdgeMesh
+ (
+ const IOobject& io,
+ const pointField& pts,
+ const edgeList& eds,
+ label concaveStart,
+ label mixedStart,
+ label nonFeatureStart,
+ label internalStart,
+ label flatStart,
+ label openStart,
+ label multipleStart,
+ const vectorField& normals,
+ const vectorField& edgeDirections,
+ const labelListList& edgeNormals,
+ const labelListList& featurePointNormals,
+ const labelList& regionEdges
+ );
+
+
+ //- Destructor
+ ~featureEdgeMesh();
+
+
+ // Member Functions
+
+ // Find
+
+ //- Find nearest surface edge for the sample point.
+ void nearestFeatureEdge
+ (
+ const point& sample,
+ scalar searchDistSqr,
+ pointIndexHit& info
+ ) const;
+
+ //- Find nearest surface edge for each sample point.
+ void nearestFeatureEdge
+ (
+ const pointField& samples,
+ const scalarField& searchDistSqr,
+ List<pointIndexHit>& info
+ ) const;
+
+ //- Find the nearest point on each type of feature edge
+ void nearestFeatureEdgeByType
+ (
+ const point& sample,
+ const scalarField& searchDistSqr,
+ List<pointIndexHit>& info
+ ) const;
+
+ // Access
+
+ //- Return the index of the start of the convex feature points
+ inline label convexStart() const;
+
+ //- Return the index of the start of the concave feature points
+ inline label concaveStart() const;
+
+ //- Return the index of the start of the mixed type feature points
+ inline label mixedStart() const;
+
+ //- Return the index of the start of the non-feature points
+ inline label nonFeatureStart() const;
+
+ //- Return the index of the start of the external feature edges
+ inline label externalStart() const;
+
+ //- Return the index of the start of the internal feature edges
+ inline label internalStart() const;
+
+ //- Return the index of the start of the flat feature edges
+ inline label flatStart() const;
+
+ //- Return the index of the start of the open feature edges
+ inline label openStart() const;
+
+ //- Return the index of the start of the multiply-connected feature
+ // edges
+ inline label multipleStart() const;
+
+ //- Return whether or not the point index is a feature point
+ inline bool featurePoint(label ptI) const;
- //- WriteData function required for regIOobject write operation
- virtual bool writeData(Ostream&) const;
+ //- Return the normals of the surfaces adjacent to the feature edges
+ // and points
+ inline const vectorField& normals() const;
+
+ //- Return the edgeDirection vectors
+ inline const vectorField& edgeDirections() const;
+
+ //- Return the direction of edgeI, pointing away from ptI
+ inline vector edgeDirection(label edgeI, label ptI) const;
+
+ //- Return the indices of the normals that are adjacent to the
+ // feature edges
+ inline const labelListList& edgeNormals() const;
+
+ //- Return the normal vectors for a given set of normal indices
+ inline vectorField edgeNormals(const labelList& edgeNormIs) const;
+
+ //- Return the normal vectors for a given edge
+ inline vectorField edgeNormals(label edgeI) const;
+
+ //- Return the indices of the normals that are adjacent to the
+ // feature points
+ inline const labelListList& featurePointNormals() const;
+
+ //- Return the normal vectors for a given feature point
+ inline vectorField featurePointNormals(label ptI) const;
+
+ //- Return the feature edges which are on the boundary between
+ // regions
+ inline const labelList& regionEdges() const;
+
+ //- Return the pointStatus of a specified point
+ inline pointStatus getPointStatus(label ptI) const;
+
+ //- Return the edgeStatus of a specified edge
+ inline edgeStatus getEdgeStatus(label edgeI) const;
+
+ //- Demand driven construction of octree for boundary edges
+ const indexedOctree<treeDataEdge>& edgeTree() const;
+
+ //- Demand driven construction of octree for boundary edges by type
+ const PtrList<indexedOctree<treeDataEdge> >&
+ edgeTreesByType() const;
+
+
+ // Write
+
+ //- Write all components of the featureEdgeMesh as obj files
+ void writeObj(const fileName& prefix) const;
+
+ //- Give precedence to the regIOobject write
+ using regIOobject::write;
+
+ //- WriteData function required for regIOobject write operation
+ virtual bool writeData(Ostream&) const;
};
@@ -92,6 +362,10 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#include "featureEdgeMeshI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
#endif
// ************************************************************************* //
--
GitLab