diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index c2baec519e8cbd589e3fa986be039f9c923ae829..bf44c1cd77d86b75322959ff61a62074321ba43c 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -74,30 +74,6 @@ scalarField calcCurvature(const triSurface& surf)
// The rest of this function adapted from
// CGAL-3.7/examples/Jet_fitting_3/Mesh_estimation.cpp
- // Licensed under CGAL-3.7/LICENSE.FREE_USE
-
- // Copyright (c) 1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007
- // Utrecht University (The Netherlands), ETH Zurich (Switzerland), Freie
- // Universitaet Berlin (Germany), INRIA Sophia-Antipolis (France),
- // Martin-Luther-University Halle-Wittenberg (Germany), Max-Planck-Institute
- // Saarbruecken (Germany), RISC Linz (Austria), and Tel-Aviv University
- // (Israel). All rights reserved.
-
- // Permission is hereby granted, free of charge, to any person obtaining a
- // copy of this software and associated documentation files (the
- // "Software"), to deal in the Software without restriction, including
- // without limitation the rights to use, copy, modify, merge, publish,
- // distribute, sublicense, and/or sell copies of the Software, and to permit
- // persons to whom the Software is furnished to do so, subject to the
- // following conditions:
-
- // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
- // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
- // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
- // IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
- // CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT
- // OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
- // THE USE OR OTHER DEALINGS IN THE SOFTWARE.
//Vertex property map, with std::map
typedef std::map<Vertex*, int> Vertex2int_map_type;
@@ -503,358 +479,344 @@ int main(int argc, char *argv[])
"extract and write surface features to file"
);
argList::noParallel();
- argList::validArgs.append("surface");
- argList::validArgs.append("output set");
argList::addOption
(
- "includedAngle",
- "degrees",
- "construct feature set from included angle [0..180]"
+ "dict",
+ "word",
+ "name of dictionary to provide feature extraction information"
);
- argList::addOption
- (
- "set",
- "name",
- "use existing feature set from file"
- );
- argList::addOption
- (
- "minLen",
- "scalar",
- "remove features shorter than the specified cumulative length"
- );
- argList::addOption
- (
- "minElem",
- "int",
- "remove features with fewer than the specified number of edges"
- );
- argList::addOption
- (
- "subsetBox",
- "((x0 y0 z0)(x1 y1 z1))",
- "remove edges outside specified bounding box"
- );
- argList::addOption
- (
- "deleteBox",
- "((x0 y0 z0)(x1 y1 z1))",
- "remove edges within specified bounding box"
- );
- argList::addBoolOption
- (
- "writeObj",
- "write extendedFeatureEdgeMesh obj files"
- );
- argList::addBoolOption
- (
- "writeVTK",
- "write extendedFeatureEdgeMesh vtk files"
- );
- argList::addBoolOption
- (
- "closeness",
- "span to look for surface closeness"
- );
- argList::addOption
- (
- "featureProximity",
- "scalar",
- "distance to look for close features"
- );
- argList::addBoolOption
- (
- "writeVTK",
- "write surface property VTK files"
- );
- argList::addBoolOption
- (
- "manifoldEdgesOnly",
- "remove any non-manifold (open or more than two connected faces) edges"
- );
- argList::addOption
+
+# include "setRootCase.H"
+# include "createTime.H"
+
+ word dictName
(
- "plane",
- "(nx ny nz)(z0 y0 z0)",
- "use a plane to create feature edges for 2D meshing"
+ args.optionLookupOrDefault<word>("dict", "surfaceFeatureExtractDict")
);
-# ifdef ENABLE_CURVATURE
- argList::addBoolOption
+ Info<< "Reading " << dictName << nl << endl;
+
+ IOdictionary dict
(
- "curvature",
- "calculate curvature and closeness fields"
+ IOobject
+ (
+ dictName,
+ runTime.system(),
+ runTime,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE,
+ false
+ )
);
-# endif
+ forAllConstIter(dictionary, dict, iter)
+ {
+ dictionary surfaceDict = dict.subDict(iter().keyword());
-# include "setRootCase.H"
-# include "createTime.H"
+ const fileName surfFileName = iter().keyword();
+ const fileName sFeatFileName = surfFileName.lessExt().name();
- bool writeVTK = args.optionFound("writeVTK");
+ Info<< "Surface : " << surfFileName << nl << endl;
- bool writeObj = args.optionFound("writeObj");
+ const Switch writeVTK =
+ surfaceDict.lookupOrAddDefault<Switch>("writeVTK", "off");
+ const Switch writeObj =
+ surfaceDict.lookupOrAddDefault<Switch>("writeObj", "off");
+ const Switch writeFeatureEdgeMesh =
+ surfaceDict.lookupOrAddDefault<Switch>
+ (
+ "writeFeatureEdgeMesh",
+ "off"
+ );
- bool curvature = args.optionFound("curvature");
+ const Switch curvature =
+ surfaceDict.lookupOrAddDefault<Switch>("curvature", "off");
+ const Switch featureProximity =
+ surfaceDict.lookupOrAddDefault<Switch>("featureProximity", "off");
+ const Switch closeness =
+ surfaceDict.lookupOrAddDefault<Switch>("closeness", "off");
- if (curvature && env("FOAM_SIGFPE"))
- {
- WarningIn(args.executable())
- << "Detected floating point exception trapping (FOAM_SIGFPE)."
- << " This might give" << nl
- << " problems when calculating curvature on straight angles"
- << " (infinite curvature)" << nl
- << " Switch it off in case of problems." << endl;
- }
+ const word extractionMethod = surfaceDict.lookup("extractionMethod");
- Info<< "Feature line extraction is only valid on closed manifold surfaces."
- << endl;
+#ifndef ENABLE_CURVATURE
+ if (curvature)
+ {
+ WarningIn(args.executable())
+ << "Curvature calculation has been requested but "
+ << args.executable() << " has not " << nl
+ << " been compiled with CGAL. "
+ << "Skipping the curvature calculation." << endl;
+ }
+#else
+ if (curvature && env("FOAM_SIGFPE"))
+ {
+ WarningIn(args.executable())
+ << "Detected floating point exception trapping (FOAM_SIGFPE)."
+ << " This might give" << nl
+ << " problems when calculating curvature on straight angles"
+ << " (infinite curvature)" << nl
+ << " Switch it off in case of problems." << endl;
+ }
+#endif
- const fileName surfFileName = args[1];
- const fileName outFileName = args[2];
- Info<< "Surface : " << surfFileName << nl
- << "Output feature set : " << outFileName << nl
- << endl;
+ Info<< nl << "Feature line extraction is only valid on closed manifold "
+ << "surfaces." << endl;
- fileName sFeatFileName = surfFileName.lessExt().name();
+ // Read
+ // ~~~~
+ triSurface surf("constant/triSurface/" + surfFileName);
- // Read
- // ~~~~
+ Info<< "Statistics:" << endl;
+ surf.writeStats(Info);
+ Info<< endl;
- triSurface surf(surfFileName);
+ faceList faces(surf.size());
- Info<< "Statistics:" << endl;
- surf.writeStats(Info);
- Info<< endl;
+ forAll(surf, fI)
+ {
+ faces[fI] = surf[fI].triFaceFace();
+ }
- faceList faces(surf.size());
- forAll(surf, fI)
- {
- faces[fI] = surf[fI].triFaceFace();
- }
+ // Either construct features from surface & featureAngle or read set.
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- // Either construct features from surface&featureangle or read set.
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ surfaceFeatures set(surf);
- surfaceFeatures set(surf);
+ if (extractionMethod == "extractFromFile")
+ {
+ const fileName featureEdgeFile =
+ surfaceDict.subDict("extractFromFile").lookup
+ (
+ "featureEdgeFile"
+ );
- if (args.optionFound("set"))
- {
- const fileName setName = args["set"];
+ edgeMesh eMesh(featureEdgeFile);
- Info<< "Reading existing feature set from file " << setName << endl;
+ // Sometimes duplicate edges are present. Remove them.
+ eMesh.mergeEdges();
- set = surfaceFeatures(surf, setName);
- }
- else if (args.optionFound("includedAngle"))
- {
- const scalar includedAngle = args.optionRead<scalar>("includedAngle");
+ Info<< nl << "Reading existing feature edges from file "
+ << featureEdgeFile << endl;
+
+ set = surfaceFeatures(surf, eMesh.points(), eMesh.edges());
+ }
+ else if (extractionMethod == "extractFromSurface")
+ {
+ const scalar includedAngle =
+ readScalar
+ (
+ surfaceDict.subDict("extractFromSurface").lookup
+ (
+ "includedAngle"
+ )
+ );
- Info<< "Constructing feature set from included angle " << includedAngle
+ Info<< nl << "Constructing feature set from included angle "
+ << includedAngle << endl;
+
+ set = surfaceFeatures(surf, includedAngle);
+ }
+ else
+ {
+ FatalErrorIn(args.executable())
+ << "No initial feature set. Provide either one"
+ << " of extractFromFile (to read existing set)" << nl
+ << " or extractFromSurface (to construct new set from angle)"
+ << exit(FatalError);
+ }
+
+ Info<< nl
+ << "Initial feature set:" << nl
+ << " feature points : " << set.featurePoints().size() << nl
+ << " feature edges : " << set.featureEdges().size() << nl
+ << " of which" << nl
+ << " region edges : " << set.nRegionEdges() << nl
+ << " external edges : " << set.nExternalEdges() << nl
+ << " internal edges : " << set.nInternalEdges() << nl
<< endl;
- set = surfaceFeatures(surf, includedAngle);
- // Info<< nl << "Writing initial features" << endl;
- // set.write("initial.fSet");
- // set.writeObj("initial");
- }
- else
- {
- FatalErrorIn(args.executable())
- << "No initial feature set. Provide either one"
- << " of -set (to read existing set)" << nl
- << " or -includedAngle (to new set construct from angle)"
- << exit(FatalError);
- }
+ // Trim set
+ // ~~~~~~~~
- Info<< nl
- << "Initial feature set:" << nl
- << " feature points : " << set.featurePoints().size() << nl
- << " feature edges : " << set.featureEdges().size() << nl
- << " of which" << nl
- << " region edges : " << set.nRegionEdges() << nl
- << " external edges : " << set.nExternalEdges() << nl
- << " internal edges : " << set.nInternalEdges() << nl
- << endl;
+ if (surfaceDict.isDict("trimFeatures"))
+ {
+ dictionary trimDict = surfaceDict.subDict("trimFeatures");
+ scalar minLen =
+ trimDict.lookupOrAddDefault<scalar>("minLen", -GREAT);
- // Trim set
- // ~~~~~~~~
+ label minElem = trimDict.lookupOrAddDefault<label>("minElem", 0);
- scalar minLen = -GREAT;
- if (args.optionReadIfPresent("minLen", minLen))
- {
- Info<< "Removing features of length < " << minLen << endl;
- }
+ // Trim away small groups of features
+ if (minElem > 0 || minLen > 0)
+ {
+ Info<< "Removing features of length < "
+ << minLen << endl;
+ Info<< "Removing features with number of edges < "
+ << minElem << endl;
- label minElem = 0;
- if (args.optionReadIfPresent("minElem", minElem))
- {
- Info<< "Removing features with number of edges < " << minElem << endl;
- }
+ set.trimFeatures(minLen, minElem);
+ }
+ }
- // Trim away small groups of features
- if (minElem > 0 || minLen > 0)
- {
- set.trimFeatures(minLen, minElem);
- Info<< endl << "Removed small features" << endl;
- }
+ // Subset
+ // ~~~~~~
- // Subset
- // ~~~~~~
+ // Convert to marked edges, points
+ List<surfaceFeatures::edgeStatus> edgeStat(set.toStatus());
- // Convert to marked edges, points
- List<surfaceFeatures::edgeStatus> edgeStat(set.toStatus());
+ if (surfaceDict.isDict("subsetFeatures"))
+ {
+ dictionary subsetDict = surfaceDict.subDict("subsetFeatures");
- if (args.optionFound("subsetBox"))
- {
- treeBoundBox bb
- (
- args.optionLookup("subsetBox")()
- );
+ if (subsetDict.found("insideBox"))
+ {
+ treeBoundBox bb(subsetDict.lookup("insideBox")());
- Info<< "Removing all edges outside bb " << bb << endl;
- dumpBox(bb, "subsetBox.obj");
+ Info<< "Removing all edges outside bb " << bb << endl;
+ dumpBox(bb, "subsetBox.obj");
- deleteBox
- (
- surf,
- bb,
- false,
- edgeStat
- );
- }
- else if (args.optionFound("deleteBox"))
- {
- treeBoundBox bb
- (
- args.optionLookup("deleteBox")()
- );
+ deleteBox
+ (
+ surf,
+ bb,
+ false,
+ edgeStat
+ );
+ }
+ else if (subsetDict.found("outsideBox"))
+ {
+ treeBoundBox bb(subsetDict.lookup("outsideBox")());
- Info<< "Removing all edges inside bb " << bb << endl;
- dumpBox(bb, "deleteBox.obj");
+ Info<< "Removing all edges inside bb " << bb << endl;
+ dumpBox(bb, "deleteBox.obj");
- deleteBox
- (
- surf,
- bb,
- true,
- edgeStat
- );
- }
+ deleteBox
+ (
+ surf,
+ bb,
+ true,
+ edgeStat
+ );
+ }
- if (args.optionFound("manifoldEdgesOnly"))
- {
- Info<< "Removing all non-manifold edges" << endl;
+ const Switch manifoldEdges =
+ subsetDict.lookupOrAddDefault<Switch>("manifoldEdges", "no");
- forAll(edgeStat, edgeI)
- {
- if (surf.edgeFaces()[edgeI].size() != 2)
+ if (manifoldEdges)
{
- edgeStat[edgeI] = surfaceFeatures::NONE;
+ Info<< "Removing all non-manifold edges" << endl;
+
+ forAll(edgeStat, edgeI)
+ {
+ if (surf.edgeFaces()[edgeI].size() != 2)
+ {
+ edgeStat[edgeI] = surfaceFeatures::NONE;
+ }
+ }
}
- }
- }
- if (args.optionFound("plane"))
- {
- plane cutPlane(args.optionLookup("plane")());
+ if (subsetDict.found("plane"))
+ {
+ plane cutPlane(subsetDict.lookup("plane")());
- deleteEdges
- (
- surf,
- cutPlane,
- edgeStat
- );
+ deleteEdges
+ (
+ surf,
+ cutPlane,
+ edgeStat
+ );
- Info<< "Only edges that intersect the plane with normal "
- << cutPlane.normal() << " and base point " << cutPlane.refPoint()
- << " will be included as feature edges."<< endl;
- }
+ Info<< "Only edges that intersect the plane with normal "
+ << cutPlane.normal()
+ << " and base point " << cutPlane.refPoint()
+ << " will be included as feature edges."<< endl;
+ }
+ }
- surfaceFeatures newSet(surf);
- newSet.setFromStatus(edgeStat);
+ surfaceFeatures newSet(surf);
+ newSet.setFromStatus(edgeStat);
- //Info<< endl << "Writing trimmed features to " << outFileName << endl;
- //newSet.write(outFileName);
+ if (writeObj)
+ {
+ newSet.writeObj("final");
+ }
- // Info<< endl << "Writing edge objs." << endl;
- // newSet.writeObj("final");
+ Info<< nl
+ << "Final feature set after trimming and subsetting:" << nl
+ << " feature points : " << newSet.featurePoints().size() << nl
+ << " feature edges : " << newSet.featureEdges().size() << nl
+ << " of which" << nl
+ << " region edges : " << newSet.nRegionEdges() << nl
+ << " external edges : " << newSet.nExternalEdges() << nl
+ << " internal edges : " << newSet.nInternalEdges() << nl
+ << endl;
- Info<< nl
- << "Final feature set:" << nl
- << " feature points : " << newSet.featurePoints().size() << nl
- << " feature edges : " << newSet.featureEdges().size() << nl
- << " of which" << nl
- << " region edges : " << newSet.nRegionEdges() << nl
- << " external edges : " << newSet.nExternalEdges() << nl
- << " internal edges : " << newSet.nInternalEdges() << nl
- << endl;
+ // Extracting and writing a extendedFeatureEdgeMesh
+ extendedFeatureEdgeMesh feMesh
+ (
+ newSet,
+ runTime,
+ sFeatFileName + ".extendedFeatureEdgeMesh"
+ );
- // Extracting and writing a extendedFeatureEdgeMesh
- extendedFeatureEdgeMesh feMesh
- (
- newSet,
- runTime,
- sFeatFileName + ".extendedFeatureEdgeMesh"
- );
+ Info<< nl << "Writing extendedFeatureEdgeMesh to "
+ << feMesh.objectPath() << endl;
- Info<< nl << "Writing extendedFeatureEdgeMesh to " << feMesh.objectPath()
- << endl;
+ if (writeObj)
+ {
+ feMesh.writeObj(surfFileName.lessExt().name());
+ }
- if (writeObj)
- {
- feMesh.writeObj(surfFileName.lessExt().name());
- }
+ feMesh.write();
- feMesh.write();
+ // Write a featureEdgeMesh for backwards compatibility
+ if (writeFeatureEdgeMesh)
+ {
+ featureEdgeMesh bfeMesh
+ (
+ IOobject
+ (
+ surfFileName.lessExt().name() + ".eMesh", // name
+ runTime.constant(), // instance
+ "triSurface",
+ runTime, // registry
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE,
+ false
+ ),
+ feMesh.points(),
+ feMesh.edges()
+ );
- // Write a featureEdgeMesh for backwards compatibility
- {
- featureEdgeMesh bfeMesh
+ Info<< nl << "Writing featureEdgeMesh to "
+ << bfeMesh.objectPath() << endl;
+
+ bfeMesh.regIOobject::write();
+ }
+
+ triSurfaceMesh searchSurf
(
IOobject
(
- surfFileName.lessExt().name() + ".eMesh", // name
- runTime.constant(), // instance
+ sFeatFileName + ".closeness",
+ runTime.constant(),
"triSurface",
- runTime, // registry
+ runTime,
IOobject::NO_READ,
- IOobject::AUTO_WRITE,
- false
+ IOobject::NO_WRITE
),
- feMesh.points(),
- feMesh.edges()
+ surf
);
- Info<< nl << "Writing featureEdgeMesh to "
- << bfeMesh.objectPath() << endl;
-
- bfeMesh.regIOobject::write();
- }
-
- triSurfaceMesh searchSurf
- (
- IOobject
- (
- sFeatFileName + ".closeness",
- runTime.constant(),
- "triSurface",
- runTime,
- IOobject::NO_READ,
- IOobject::NO_WRITE
- ),
- surf
- );
-
// Find close features
// // Dummy trim operation to mark features
@@ -922,72 +884,59 @@ int main(int argc, char *argv[])
// )
// );
- if (args.optionFound("closeness"))
- {
- Info<< nl << "Extracting internal and external closeness of surface."
- << endl;
-
- // Internal and external closeness
+ if (closeness)
+ {
+ Info<< nl << "Extracting internal and external closeness of "
+ << "surface." << endl;
- // Prepare start and end points for intersection tests
+ // Internal and external closeness
- const vectorField& normals = searchSurf.faceNormals();
+ // Prepare start and end points for intersection tests
- scalar span = searchSurf.bounds().mag();
+ const vectorField& normals = searchSurf.faceNormals();
- //args.optionReadIfPresent("closeness", span);
+ scalar span = searchSurf.bounds().mag();
- scalar externalAngleTolerance = 10;
- scalar externalToleranceCosAngle = Foam::cos
- (
- degToRad(180 - externalAngleTolerance)
- );
+ //args.optionReadIfPresent("closeness", span);
- scalar internalAngleTolerance = 45;
- scalar internalToleranceCosAngle = Foam::cos
- (
- degToRad(180 - internalAngleTolerance)
- );
+ scalar externalAngleTolerance = 10;
+ scalar externalToleranceCosAngle =
+ Foam::cos
+ (
+ degToRad(180 - externalAngleTolerance)
+ );
- Info<< "externalToleranceCosAngle: " << externalToleranceCosAngle << nl
- << "internalToleranceCosAngle: " << internalToleranceCosAngle
- << endl;
+ scalar internalAngleTolerance = 45;
+ scalar internalToleranceCosAngle =
+ Foam::cos
+ (
+ degToRad(180 - internalAngleTolerance)
+ );
- // Info<< "span " << span << endl;
+ Info<< "externalToleranceCosAngle: " << externalToleranceCosAngle
+ << nl
+ << "internalToleranceCosAngle: " << internalToleranceCosAngle
+ << endl;
- pointField start = searchSurf.faceCentres() - span*normals;
- pointField end = searchSurf.faceCentres() + span*normals;
- const pointField& faceCentres = searchSurf.faceCentres();
+ // Info<< "span " << span << endl;
- List<List<pointIndexHit> > allHitInfo;
+ pointField start = searchSurf.faceCentres() - span*normals;
+ pointField end = searchSurf.faceCentres() + span*normals;
+ const pointField& faceCentres = searchSurf.faceCentres();
- // Find all intersections (in order)
- searchSurf.findLineAll(start, end, allHitInfo);
+ List<List<pointIndexHit> > allHitInfo;
- scalarField internalCloseness(start.size(), GREAT);
- scalarField externalCloseness(start.size(), GREAT);
+ // Find all intersections (in order)
+ searchSurf.findLineAll(start, end, allHitInfo);
- forAll(allHitInfo, fI)
- {
- const List<pointIndexHit>& hitInfo = allHitInfo[fI];
+ scalarField internalCloseness(start.size(), GREAT);
+ scalarField externalCloseness(start.size(), GREAT);
- if (hitInfo.size() < 1)
+ forAll(allHitInfo, fI)
{
- drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
+ const List<pointIndexHit>& hitInfo = allHitInfo[fI];
- // FatalErrorIn(args.executable())
- // << "findLineAll did not hit its own face."
- // << exit(FatalError);
- }
- else if (hitInfo.size() == 1)
- {
- if (!hitInfo[0].hit())
- {
- // FatalErrorIn(args.executable())
- // << "findLineAll did not hit any face."
- // << exit(FatalError);
- }
- else if (hitInfo[0].index() != fI)
+ if (hitInfo.size() < 1)
{
drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
@@ -995,283 +944,339 @@ int main(int argc, char *argv[])
// << "findLineAll did not hit its own face."
// << exit(FatalError);
}
- }
- else
- {
- label ownHitI = -1;
-
- forAll(hitInfo, hI)
+ else if (hitInfo.size() == 1)
{
- // Find the hit on the triangle that launched the ray
-
- if (hitInfo[hI].index() == fI)
+ if (!hitInfo[0].hit())
{
- ownHitI = hI;
-
- break;
+ // FatalErrorIn(args.executable())
+ // << "findLineAll did not hit any face."
+ // << exit(FatalError);
}
- }
-
- if (ownHitI < 0)
- {
- drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
-
- // FatalErrorIn(args.executable())
- // << "findLineAll did not hit its own face."
- // << exit(FatalError);
- }
- else if (ownHitI == 0)
- {
- // There are no internal hits, the first hit is the closest
- // external hit
-
- if
- (
- (normals[fI] & normals[hitInfo[ownHitI + 1].index()])
- < externalToleranceCosAngle
- )
+ else if (hitInfo[0].index() != fI)
{
- externalCloseness[fI] = mag
+ drawHitProblem
(
- faceCentres[fI] - hitInfo[ownHitI + 1].hitPoint()
+ fI,
+ surf,
+ start,
+ faceCentres,
+ end,
+ hitInfo
);
+
+ // FatalErrorIn(args.executable())
+ // << "findLineAll did not hit its own face."
+ // << exit(FatalError);
}
}
- else if (ownHitI == hitInfo.size() - 1)
+ else
{
- // There are no external hits, the last but one hit is the
- // closest internal hit
+ label ownHitI = -1;
- if
- (
- (normals[fI] & normals[hitInfo[ownHitI - 1].index()])
- < internalToleranceCosAngle
- )
+ forAll(hitInfo, hI)
+ {
+ // Find the hit on the triangle that launched the ray
+
+ if (hitInfo[hI].index() == fI)
+ {
+ ownHitI = hI;
+
+ break;
+ }
+ }
+
+ if (ownHitI < 0)
{
- internalCloseness[fI] = mag
+ drawHitProblem
(
- faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint()
+ fI,
+ surf,
+ start,
+ faceCentres,
+ end,
+ hitInfo
);
+
+ // FatalErrorIn(args.executable())
+ // << "findLineAll did not hit its own face."
+ // << exit(FatalError);
}
- }
- else
- {
- if
- (
- (normals[fI] & normals[hitInfo[ownHitI + 1].index()])
- < externalToleranceCosAngle
- )
+ else if (ownHitI == 0)
{
- externalCloseness[fI] = mag
+ // There are no internal hits, the first hit is the
+ // closest external hit
+
+ if
(
- faceCentres[fI] - hitInfo[ownHitI + 1].hitPoint()
- );
+ (
+ normals[fI]
+ & normals[hitInfo[ownHitI + 1].index()]
+ )
+ < externalToleranceCosAngle
+ )
+ {
+ externalCloseness[fI] =
+ mag
+ (
+ faceCentres[fI]
+ - hitInfo[ownHitI + 1].hitPoint()
+ );
+ }
}
+ else if (ownHitI == hitInfo.size() - 1)
+ {
+ // There are no external hits, the last but one hit is
+ // the closest internal hit
- if
- (
- (normals[fI] & normals[hitInfo[ownHitI - 1].index()])
- < internalToleranceCosAngle
- )
+ if
+ (
+ (
+ normals[fI]
+ & normals[hitInfo[ownHitI - 1].index()]
+ )
+ < internalToleranceCosAngle
+ )
+ {
+ internalCloseness[fI] =
+ mag
+ (
+ faceCentres[fI]
+ - hitInfo[ownHitI - 1].hitPoint()
+ );
+ }
+ }
+ else
{
- internalCloseness[fI] = mag
+ if
(
- faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint()
- );
+ (
+ normals[fI]
+ & normals[hitInfo[ownHitI + 1].index()]
+ )
+ < externalToleranceCosAngle
+ )
+ {
+ externalCloseness[fI] =
+ mag
+ (
+ faceCentres[fI]
+ - hitInfo[ownHitI + 1].hitPoint()
+ );
+ }
+
+ if
+ (
+ (
+ normals[fI]
+ & normals[hitInfo[ownHitI - 1].index()]
+ )
+ < internalToleranceCosAngle
+ )
+ {
+ internalCloseness[fI] =
+ mag
+ (
+ faceCentres[fI]
+ - hitInfo[ownHitI - 1].hitPoint()
+ );
+ }
}
}
}
- }
- triSurfaceScalarField internalClosenessField
- (
- IOobject
+ triSurfaceScalarField internalClosenessField
(
- sFeatFileName + ".internalCloseness",
- runTime.constant(),
- "triSurface",
- runTime,
- IOobject::NO_READ,
- IOobject::NO_WRITE
- ),
- surf,
- dimLength,
- internalCloseness
- );
+ IOobject
+ (
+ sFeatFileName + ".internalCloseness",
+ runTime.constant(),
+ "triSurface",
+ runTime,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ surf,
+ dimLength,
+ internalCloseness
+ );
- internalClosenessField.write();
+ internalClosenessField.write();
- triSurfaceScalarField externalClosenessField
- (
- IOobject
+ triSurfaceScalarField externalClosenessField
(
- sFeatFileName + ".externalCloseness",
- runTime.constant(),
- "triSurface",
- runTime,
- IOobject::NO_READ,
- IOobject::NO_WRITE
- ),
- surf,
- dimLength,
- externalCloseness
- );
+ IOobject
+ (
+ sFeatFileName + ".externalCloseness",
+ runTime.constant(),
+ "triSurface",
+ runTime,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ surf,
+ dimLength,
+ externalCloseness
+ );
- externalClosenessField.write();
+ externalClosenessField.write();
- if (writeVTK)
- {
- vtkSurfaceWriter().write
- (
- runTime.constant()/"triSurface", // outputDir
- sFeatFileName, // surfaceName
- surf.points(),
- faces,
- "internalCloseness", // fieldName
- internalCloseness,
- false, // isNodeValues
- true // verbose
- );
+ if (writeVTK)
+ {
+ vtkSurfaceWriter().write
+ (
+ runTime.constant()/"triSurface", // outputDir
+ sFeatFileName, // surfaceName
+ surf.points(),
+ faces,
+ "internalCloseness", // fieldName
+ internalCloseness,
+ false, // isNodeValues
+ true // verbose
+ );
- vtkSurfaceWriter().write
- (
- runTime.constant()/"triSurface", // outputDir
- sFeatFileName, // surfaceName
- surf.points(),
- faces,
- "externalCloseness", // fieldName
- externalCloseness,
- false, // isNodeValues
- true // verbose
- );
+ vtkSurfaceWriter().write
+ (
+ runTime.constant()/"triSurface", // outputDir
+ sFeatFileName, // surfaceName
+ surf.points(),
+ faces,
+ "externalCloseness", // fieldName
+ externalCloseness,
+ false, // isNodeValues
+ true // verbose
+ );
+ }
}
- }
#ifdef ENABLE_CURVATURE
- if (args.optionFound("curvature"))
- {
- Info<< nl << "Extracting curvature of surface at the points." << endl;
+ if (curvature)
+ {
+ Info<< nl << "Extracting curvature of surface at the points."
+ << endl;
- scalarField k = calcCurvature(surf);
+ scalarField k = calcCurvature(surf);
// Modify the curvature values on feature edges and points to be zero.
- // forAll(newSet.featureEdges(), fEI)
- // {
- // const edge& e = surf.edges()[newSet.featureEdges()[fEI]];
- //
- // k[surf.meshPoints()[e.start()]] = 0.0;
- // k[surf.meshPoints()[e.end()]] = 0.0;
- // }
+ // forAll(newSet.featureEdges(), fEI)
+ // {
+ // const edge& e = surf.edges()[newSet.featureEdges()[fEI]];
+ //
+ // k[surf.meshPoints()[e.start()]] = 0.0;
+ // k[surf.meshPoints()[e.end()]] = 0.0;
+ // }
- triSurfacePointScalarField kField
- (
- IOobject
+ triSurfacePointScalarField kField
(
- sFeatFileName + ".curvature",
- runTime.constant(),
- "triSurface",
- runTime,
- IOobject::NO_READ,
- IOobject::NO_WRITE
- ),
- surf,
- dimLength,
- k
- );
+ IOobject
+ (
+ sFeatFileName + ".curvature",
+ runTime.constant(),
+ "triSurface",
+ runTime,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ surf,
+ dimLength,
+ k
+ );
- kField.write();
+ kField.write();
- if (writeVTK)
- {
- vtkSurfaceWriter().write
- (
- runTime.constant()/"triSurface", // outputDir
- sFeatFileName, // surfaceName
- surf.points(),
- faces,
- "curvature", // fieldName
- k,
- true, // isNodeValues
- true // verbose
- );
+ if (writeVTK)
+ {
+ vtkSurfaceWriter().write
+ (
+ runTime.constant()/"triSurface", // outputDir
+ sFeatFileName, // surfaceName
+ surf.points(),
+ faces,
+ "curvature", // fieldName
+ k,
+ true, // isNodeValues
+ true // verbose
+ );
+ }
}
- }
#endif
- if (args.optionFound("featureProximity"))
- {
- Info<< nl << "Extracting proximity of close feature points and edges "
- << "to the surface" << endl;
+ if (featureProximity)
+ {
+ Info<< nl << "Extracting proximity of close feature points and "
+ << "edges to the surface" << endl;
- const scalar searchDistance =
- args.optionRead<scalar>("featureProximity");
+ const scalar searchDistance =
+ readScalar(surfaceDict.lookup("maxFeatureProximity"));
- const scalar radiusSqr = sqr(searchDistance);
+ const scalar radiusSqr = sqr(searchDistance);
- scalarField featureProximity(surf.size(), searchDistance);
+ scalarField featureProximity(surf.size(), searchDistance);
- forAll(surf, fI)
- {
- const triPointRef& tri = surf[fI].tri(surf.points());
- const point& triCentre = tri.circumCentre();
+ forAll(surf, fI)
+ {
+ const triPointRef& tri = surf[fI].tri(surf.points());
+ const point& triCentre = tri.circumCentre();
- List<pointIndexHit> hitList;
+ List<pointIndexHit> hitList;
- feMesh.allNearestFeatureEdges(triCentre, radiusSqr, hitList);
+ feMesh.allNearestFeatureEdges(triCentre, radiusSqr, hitList);
- featureProximity[fI] =
- calcProximityOfFeatureEdges
- (
- feMesh,
- hitList,
- featureProximity[fI]
- );
+ featureProximity[fI] =
+ calcProximityOfFeatureEdges
+ (
+ feMesh,
+ hitList,
+ featureProximity[fI]
+ );
- feMesh.allNearestFeaturePoints(triCentre, radiusSqr, hitList);
+ feMesh.allNearestFeaturePoints(triCentre, radiusSqr, hitList);
- featureProximity[fI] =
- calcProximityOfFeaturePoints
- (
- hitList,
- featureProximity[fI]
- );
- }
+ featureProximity[fI] =
+ calcProximityOfFeaturePoints
+ (
+ hitList,
+ featureProximity[fI]
+ );
+ }
- triSurfaceScalarField featureProximityField
- (
- IOobject
+ triSurfaceScalarField featureProximityField
(
- sFeatFileName + ".featureProximity",
- runTime.constant(),
- "triSurface",
- runTime,
- IOobject::NO_READ,
- IOobject::NO_WRITE
- ),
- surf,
- dimLength,
- featureProximity
- );
+ IOobject
+ (
+ sFeatFileName + ".featureProximity",
+ runTime.constant(),
+ "triSurface",
+ runTime,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ surf,
+ dimLength,
+ featureProximity
+ );
- featureProximityField.write();
+ featureProximityField.write();
- if (writeVTK)
- {
- vtkSurfaceWriter().write
- (
- runTime.constant()/"triSurface", // outputDir
- sFeatFileName, // surfaceName
- surf.points(),
- faces,
- "featureProximity", // fieldName
- featureProximity,
- false, // isNodeValues
- true // verbose
- );
+ if (writeVTK)
+ {
+ vtkSurfaceWriter().write
+ (
+ runTime.constant()/"triSurface", // outputDir
+ sFeatFileName, // surfaceName
+ surf.points(),
+ faces,
+ "featureProximity", // fieldName
+ featureProximity,
+ false, // isNodeValues
+ true // verbose
+ );
+ }
}
+
+ Info<< endl;
}
Info<< "End\n" << endl;
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict
new file mode 100644
index 0000000000000000000000000000000000000000..3621d64a293ffa8a1a2daf3be441818036555e92
--- /dev/null
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict
@@ -0,0 +1,86 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ object surfaceFeatureExtractDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+surface1.stl
+{
+ // extractFromFile || extractFromSurface
+ extractionMethod extractFromFile;
+
+ extractFromFile
+ {
+ // Load from an existing feature edge file
+ featureEdgeFile "constant/triSurface/featureEdges.nas";
+ }
+
+ trimFeatures
+ {
+ // Remove features with fewer than the specified number of edges
+ minElem 0;
+
+ // Remove features shorter than the specified cumulative length
+ minLen 0.0;
+ }
+
+ subsetFeatures
+ {
+ // Use a plane to select feature edges
+ // (normal)(basePoint)
+ plane (1 0 0)(0 0 0);
+
+ // Select feature edges using a box
+ // (minPt)(maxPt)
+ insideBox (0 0 0)(1 1 1);
+ outsideBox (0 0 0)(1 1 1);
+
+ // Remove any non-manifold (open or > 2 connected faces) edges
+ manifoldEdges no;
+ }
+
+ // Output the curvature of the surface
+ curvature no;
+
+ // Output the proximity of feature points and edges to each other
+ featureProximity no;
+ // The maximum search distance to use when looking for other feature
+ // points and edges
+ maxFeatureProximity 1;
+
+ // Out put the closeness of surface elements to other surface elements.
+ closeness no;
+
+ // Write options
+ writeVTK no;
+ writeObj yes;
+ writeFeatureEdgeMesh no;
+}
+
+
+surface2.nas
+{
+ extractionMethod extractFromSurface;
+
+ extractFromSurface
+ {
+ // Mark edges whose adjacent surface normals are at an angle less
+ // than includedAngle as features
+ // - 0 : selects no edges
+ // - 180: selects all edges
+ includedAngle 120;
+ }
+}
+
+
+// ************************************************************************* //