diff --git a/applications/utilities/surface/surfaceFeatureExtract/Make/files b/applications/utilities/surface/surfaceFeatureExtract/Make/files index f8ead594d..44b4b18fc 100644 --- a/applications/utilities/surface/surfaceFeatureExtract/Make/files +++ b/applications/utilities/surface/surfaceFeatureExtract/Make/files @@ -1,5 +1,3 @@ -surfaceExtractCloseness.C -surfaceExtractPointCloseness.C surfaceFeatureExtractUtilities.C surfaceFeatureExtract.C diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceExtractPointCloseness.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceExtractPointCloseness.C deleted file mode 100644 index 14001e662..000000000 --- a/applications/utilities/surface/surfaceFeatureExtract/surfaceExtractPointCloseness.C +++ /dev/null @@ -1,150 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2018 OpenFOAM Foundation - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM. If not, see . - -\*---------------------------------------------------------------------------*/ - -#include "surfaceFeatureExtract.H" -#include "Time.H" -#include "vtkSurfaceWriter.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -Foam::Pair> -Foam::extractPointCloseness -( - const triSurfaceMesh& surf -) -{ - const Time& runTime = surf.objectRegistry::time(); - - // Prepare start and end points for intersection tests - - const pointField& points = surf.points(); - const labelList& meshPoints = surf.meshPoints(); - const pointField& faceCentres = surf.faceCentres(); - const vectorField& normals = surf.faceNormals(); - const labelListList& pointFaces = surf.pointFaces(); - - const scalar span = surf.bounds().mag(); - - label nPointFaces = 0; - forAll(pointFaces, pfi) - { - nPointFaces += pointFaces[pfi].size(); - } - - pointField facePoints(nPointFaces); - pointField start(nPointFaces); - pointField end(nPointFaces); - - label i = 0; - forAll(points, pi) - { - forAll(pointFaces[pi], pfi) - { - const label fi = pointFaces[pi][pfi]; - - facePoints[i] = (0.9*points[meshPoints[pi]] + 0.1*faceCentres[fi]); - const vector& n = normals[fi]; - - start[i] = facePoints[i] - span*n; - end[i] = facePoints[i] + span*n; - - i++; - } - } - - List allHitinfo; - - // Find all intersections (in order) - surf.findLineAll(start, end, allHitinfo); - - scalarField internalCloseness(points.size(), great); - scalarField externalCloseness(points.size(), great); - - i = 0; - forAll(points, pi) - { - forAll(pointFaces[pi], pfi) - { - const label fi = pointFaces[pi][pfi]; - const pointIndexHitList& hitInfo = allHitinfo[i]; - - processHit - ( - internalCloseness[pi], - externalCloseness[pi], - fi, - surf, - start[i], - facePoints[i], - end[i], - normals[fi], - normals, - hitInfo - ); - - i++; - } - } - - return Pair> - ( - tmp - ( - new triSurfacePointScalarField - ( - IOobject - ( - surf.objectRegistry::name() + ".internalPointCloseness", - runTime.constant(), - "triSurface", - runTime - ), - surf, - dimLength, - internalCloseness - ) - ), - - tmp - ( - new triSurfacePointScalarField - ( - IOobject - ( - surf.objectRegistry::name() + ".externalPointCloseness", - runTime.constant(), - "triSurface", - runTime - ), - surf, - dimLength, - externalCloseness - ) - ) - ); -} - - -// ************************************************************************* // diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C index bfe9c6e5e..17f79f3be 100644 --- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C +++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C @@ -399,11 +399,6 @@ int main(int argc, char *argv[]) Info<< nl << "Extracting internal and external closeness of " << "surface." << endl; - Info<< "externalToleranceCosAngle: " << externalToleranceCosAngle - << nl - << "internalToleranceCosAngle: " << internalToleranceCosAngle - << endl; - // Searchable triSurface const triSurfaceMesh searchSurf ( @@ -420,7 +415,7 @@ int main(int argc, char *argv[]) { Pair> closenessFields ( - extractCloseness(searchSurf) + searchSurf.extractCloseness() ); closenessFields.first()->write(); @@ -459,7 +454,7 @@ int main(int argc, char *argv[]) { Pair> closenessFields ( - extractPointCloseness(searchSurf) + searchSurf.extractPointCloseness() ); closenessFields.first()->write(); diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.H b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.H index 7a1e5c7d9..1ce1a4eab 100644 --- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.H +++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.H @@ -39,12 +39,6 @@ Description namespace Foam { - extern const scalar internalAngleTolerance; - extern const scalar internalToleranceCosAngle; - - extern const scalar externalAngleTolerance; - extern const scalar externalToleranceCosAngle; - //- Deletes all edges inside/outside bounding box from set. void deleteBox ( @@ -62,18 +56,6 @@ namespace Foam List& edgeStat ); - //- Divide into multiple normal bins - // - return REGION if != 2 normals - // - return REGION if 2 normals that make feature angle - // - otherwise return NONE and set normals,bins - surfaceFeatures::edgeStatus checkNonManifoldEdge - ( - const triSurface& surf, - const scalar tol, - const scalar includedAngle, - const label edgei - ); - void deleteNonManifoldEdges ( const triSurface& surf, @@ -83,41 +65,6 @@ namespace Foam ); void writeStats(const extendedFeatureEdgeMesh& fem, Ostream& os); - - - void drawHitProblem - ( - const label fi, - const triSurface& surf, - const point& start, - const point& p, - const point& end, - const pointIndexHitList& hitInfo - ); - - void processHit - ( - scalar& internalCloseness, - scalar& externalCloseness, - const label fi, - const triSurface& surf, - const point& start, - const point& p, - const point& end, - const vector& normal, - const vectorField& normals, - const pointIndexHitList& hitInfo - ); - - Pair> extractCloseness - ( - const triSurfaceMesh& surf - ); - - Pair> extractPointCloseness - ( - const triSurfaceMesh& surf - ); } diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractUtilities.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractUtilities.C index 3d0673e71..5afe34da7 100644 --- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractUtilities.C +++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractUtilities.C @@ -87,166 +87,6 @@ void Foam::deleteEdges } -Foam::surfaceFeatures::edgeStatus Foam::checkNonManifoldEdge -( - const triSurface& surf, - const scalar tol, - const scalar includedAngle, - const label edgei -) -{ - const edge& e = surf.edges()[edgei]; - const labelList& eFaces = surf.edgeFaces()[edgei]; - - // Bin according to normal - - DynamicList normals(2); - DynamicList bins(2); - - forAll(eFaces, eFacei) - { - const Foam::vector& n = surf.faceNormals()[eFaces[eFacei]]; - - // Find the normal in normals - label index = -1; - forAll(normals, normalI) - { - if (mag(n&normals[normalI]) > (1-tol)) - { - index = normalI; - break; - } - } - - if (index != -1) - { - bins[index].append(eFacei); - } - else if (normals.size() >= 2) - { - // Would be third normal. Mark as feature. - //Pout<< "** at edge:" << surf.localPoints()[e[0]] - // << surf.localPoints()[e[1]] - // << " have normals:" << normals - // << " and " << n << endl; - return surfaceFeatures::REGION; - } - else - { - normals.append(n); - bins.append(labelList(1, eFacei)); - } - } - - - // Check resulting number of bins - if (bins.size() == 1) - { - // Note: should check here whether they are two sets of faces - // that are planar or indeed 4 faces al coming together at an edge. - //Pout<< "** at edge:" - // << surf.localPoints()[e[0]] - // << surf.localPoints()[e[1]] - // << " have single normal:" << normals[0] - // << endl; - return surfaceFeatures::NONE; - } - else - { - // Two bins. Check if normals make an angle - - //Pout<< "** at edge:" - // << surf.localPoints()[e[0]] - // << surf.localPoints()[e[1]] << nl - // << " normals:" << normals << nl - // << " bins :" << bins << nl - // << endl; - - if (includedAngle >= 0) - { - scalar minCos = Foam::cos(degToRad(180.0 - includedAngle)); - - forAll(eFaces, i) - { - const Foam::vector& ni = surf.faceNormals()[eFaces[i]]; - for (label j=i+1; j 0) - { - regionAndNormal[i] = t.region()+1; - } - else if (dir == 0) - { - FatalErrorInFunction - << exit(FatalError); - } - else - { - regionAndNormal[i] = -(t.region()+1); - } - } - - // 2. check against bin1 - const labelList& bin1 = bins[1]; - labelList regionAndNormal1(bin1.size()); - forAll(bin1, i) - { - const labelledTri& t = surf.localFaces()[eFaces[bin1[i]]]; - int dir = t.edgeDirection(e); - - label myRegionAndNormal; - if (dir > 0) - { - myRegionAndNormal = t.region()+1; - } - else - { - myRegionAndNormal = -(t.region()+1); - } - - regionAndNormal1[i] = myRegionAndNormal; - - label index = findIndex(regionAndNormal, -myRegionAndNormal); - if (index == -1) - { - // Not found. - //Pout<< "cannot find region " << myRegionAndNormal - // << " in regions " << regionAndNormal << endl; - - return surfaceFeatures::REGION; - } - } - - return surfaceFeatures::NONE; - } -} void Foam::deleteNonManifoldEdges diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files index e2c9c1d2c..4d450204a 100644 --- a/src/meshTools/Make/files +++ b/src/meshTools/Make/files @@ -114,6 +114,7 @@ searchableSurfaces/searchableSurfaces/searchableSurfaces.C searchableSurfaces/searchableSurfacesQueries/searchableSurfacesQueries.C searchableSurfaces/searchableSurfaceWithGaps/searchableSurfaceWithGaps.C searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C +searchableSurfaces/triSurfaceMesh/extractCloseness.C searchableSurfaces/closedTriSurfaceMesh/closedTriSurfaceMesh.C searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.C diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceExtractCloseness.C b/src/meshTools/searchableSurfaces/triSurfaceMesh/extractCloseness.C similarity index 53% rename from applications/utilities/surface/surfaceFeatureExtract/surfaceExtractCloseness.C rename to src/meshTools/searchableSurfaces/triSurfaceMesh/extractCloseness.C index d14a1e5fc..b5d45950a 100644 --- a/applications/utilities/surface/surfaceFeatureExtract/surfaceExtractCloseness.C +++ b/src/meshTools/searchableSurfaces/triSurfaceMesh/extractCloseness.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -23,36 +23,26 @@ License \*---------------------------------------------------------------------------*/ -#include "surfaceFeatureExtract.H" -#include "Time.H" #include "triSurfaceMesh.H" -#include "vtkSurfaceWriter.H" +#include "triSurfaceFields.H" #include "meshTools.H" +#include "Time.H" +#include "unitConversion.H" -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -const Foam::scalar Foam::internalAngleTolerance(80); -const Foam::scalar Foam::internalToleranceCosAngle -( - cos(degToRad(180 - internalAngleTolerance)) -); - -const Foam::scalar Foam::externalAngleTolerance(10); -const Foam::scalar Foam::externalToleranceCosAngle -( - cos(degToRad(180 - externalAngleTolerance)) -); - -void Foam::drawHitProblem +void Foam::triSurfaceMesh::drawHitProblem ( const label fi, - const triSurface& surf, const point& start, const point& p, const point& end, const pointIndexHitList& hitInfo -) +) const { + const List& tris = *this; + const pointField& points = this->points(); + Info<< nl << "# findLineAll did not hit its own face." << nl << "# fi " << fi << nl << "# start " << start @@ -67,19 +57,19 @@ void Foam::drawHitProblem Info<< "l 1 2 3" << endl; - meshTools::writeOBJ(Info, surf.points()[surf[fi][0]]); - meshTools::writeOBJ(Info, surf.points()[surf[fi][1]]); - meshTools::writeOBJ(Info, surf.points()[surf[fi][2]]); + meshTools::writeOBJ(Info, points[tris[fi][0]]); + meshTools::writeOBJ(Info, points[tris[fi][1]]); + meshTools::writeOBJ(Info, points[tris[fi][2]]); Info<< "f 4 5 6" << endl; forAll(hitInfo, hi) { - label hFI = hitInfo[hi].index(); + label hfi = hitInfo[hi].index(); - meshTools::writeOBJ(Info, surf.points()[surf[hFI][0]]); - meshTools::writeOBJ(Info, surf.points()[surf[hFI][1]]); - meshTools::writeOBJ(Info, surf.points()[surf[hFI][2]]); + meshTools::writeOBJ(Info, points[tris[hfi][0]]); + meshTools::writeOBJ(Info, points[tris[hfi][1]]); + meshTools::writeOBJ(Info, points[tris[hfi][2]]); Info<< "f " << 3*hi + 7 << " " @@ -90,23 +80,24 @@ void Foam::drawHitProblem } -void Foam::processHit +void Foam::triSurfaceMesh::processHit ( scalar& internalCloseness, scalar& externalCloseness, + const scalar internalToleranceCosAngle, + const scalar externalToleranceCosAngle, const label fi, - const triSurface& surf, const point& start, const point& p, const point& end, const vector& normal, const vectorField& normals, const pointIndexHitList& hitInfo -) +) const { if (hitInfo.size() < 1) { - drawHitProblem(fi, surf, start, p, end, hitInfo); + drawHitProblem(fi, start, p, end, hitInfo); } else if (hitInfo.size() == 1) { @@ -115,7 +106,7 @@ void Foam::processHit } else if (hitInfo[0].index() != fi) { - drawHitProblem(fi, surf, start, p, end, hitInfo); + drawHitProblem(fi, start, p, end, hitInfo); } } else @@ -135,7 +126,7 @@ void Foam::processHit if (ownHiti < 0) { - drawHitProblem(fi, surf, start, p, end, hitInfo); + drawHitProblem(fi, start, p, end, hitInfo); } else if (ownHiti == 0) { @@ -205,27 +196,39 @@ void Foam::processHit } -Foam::Pair> Foam::extractCloseness +Foam::Pair> +Foam::triSurfaceMesh::extractCloseness ( - const triSurfaceMesh& surf -) + const scalar internalAngleTolerance, + const scalar externalAngleTolerance +) const { - const Time& runTime = surf.objectRegistry::time(); + const scalar internalToleranceCosAngle + ( + cos(degToRad(180 - internalAngleTolerance)) + ); + + const scalar externalToleranceCosAngle + ( + cos(degToRad(180 - externalAngleTolerance)) + ); + + const Time& runTime = objectRegistry::time(); // Prepare start and end points for intersection tests - const vectorField& normals = surf.faceNormals(); + const vectorField& normals = faceNormals(); - const scalar span = surf.bounds().mag(); + const scalar span = bounds().mag(); - const pointField start(surf.faceCentres() - span*normals); - const pointField end(surf.faceCentres() + span*normals); - const pointField& faceCentres = surf.faceCentres(); + const pointField start(faceCentres() - span*normals); + const pointField end(faceCentres() + span*normals); + const pointField& faceCentres = this->faceCentres(); List> allHitinfo; // Find all intersections (in order) - surf.findLineAll(start, end, allHitinfo); + findLineAll(start, end, allHitinfo); scalarField internalCloseness(start.size(), great); scalarField externalCloseness(start.size(), great); @@ -238,8 +241,9 @@ Foam::Pair> Foam::extractCloseness ( internalCloseness[fi], externalCloseness[fi], + internalToleranceCosAngle, + externalToleranceCosAngle, fi, - surf, start[fi], faceCentres[fi], end[fi], @@ -257,12 +261,12 @@ Foam::Pair> Foam::extractCloseness ( IOobject ( - surf.objectRegistry::name() + ".internalCloseness", + objectRegistry::name() + ".internalCloseness", runTime.constant(), "triSurface", runTime ), - surf, + *this, dimLength, internalCloseness ) @@ -274,12 +278,142 @@ Foam::Pair> Foam::extractCloseness ( IOobject ( - surf.objectRegistry::name() + ".externalCloseness", + objectRegistry::name() + ".externalCloseness", runTime.constant(), "triSurface", runTime ), - surf, + *this, + dimLength, + externalCloseness + ) + ) + ); +} + + +Foam::Pair> +Foam::triSurfaceMesh::extractPointCloseness +( + const scalar internalAngleTolerance, + const scalar externalAngleTolerance +) const +{ + const scalar internalToleranceCosAngle + ( + cos(degToRad(180 - internalAngleTolerance)) + ); + + const scalar externalToleranceCosAngle + ( + cos(degToRad(180 - externalAngleTolerance)) + ); + + const Time& runTime = objectRegistry::time(); + + // Prepare start and end points for intersection tests + + const pointField& points = this->points(); + const labelList& meshPoints = this->meshPoints(); + const pointField& faceCentres = this->faceCentres(); + const vectorField& normals = this->faceNormals(); + const labelListList& pointFaces = this->pointFaces(); + + const scalar span = bounds().mag(); + + label nPointFaces = 0; + forAll(pointFaces, pfi) + { + nPointFaces += pointFaces[pfi].size(); + } + + pointField facePoints(nPointFaces); + pointField start(nPointFaces); + pointField end(nPointFaces); + + label i = 0; + forAll(points, pi) + { + forAll(pointFaces[pi], pfi) + { + const label fi = pointFaces[pi][pfi]; + + facePoints[i] = (0.9*points[meshPoints[pi]] + 0.1*faceCentres[fi]); + const vector& n = normals[fi]; + + start[i] = facePoints[i] - span*n; + end[i] = facePoints[i] + span*n; + + i++; + } + } + + List allHitinfo; + + // Find all intersections (in order) + findLineAll(start, end, allHitinfo); + + scalarField internalCloseness(points.size(), great); + scalarField externalCloseness(points.size(), great); + + i = 0; + forAll(points, pi) + { + forAll(pointFaces[pi], pfi) + { + const label fi = pointFaces[pi][pfi]; + const pointIndexHitList& hitInfo = allHitinfo[i]; + + processHit + ( + internalCloseness[pi], + externalCloseness[pi], + internalToleranceCosAngle, + externalToleranceCosAngle, + fi, + start[i], + facePoints[i], + end[i], + normals[fi], + normals, + hitInfo + ); + + i++; + } + } + + return Pair> + ( + tmp + ( + new triSurfacePointScalarField + ( + IOobject + ( + objectRegistry::name() + ".internalPointCloseness", + runTime.constant(), + "triSurface", + runTime + ), + *this, + dimLength, + internalCloseness + ) + ), + + tmp + ( + new triSurfacePointScalarField + ( + IOobject + ( + objectRegistry::name() + ".externalPointCloseness", + runTime.constant(), + "triSurface", + runTime + ), + *this, dimLength, externalCloseness ) diff --git a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C index 90d5f5f8a..72c27fbec 100644 --- a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C +++ b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C @@ -794,6 +794,37 @@ void Foam::triSurfaceMesh::getNormal } +void Foam::triSurfaceMesh::getVolumeType +( + const pointField& points, + List& volType +) const +{ + volType.setSize(points.size()); + + scalar oldTol = indexedOctree::perturbTol(); + indexedOctree::perturbTol() = tolerance(); + + forAll(points, pointi) + { + const point& pt = points[pointi]; + + if (!tree().bb().contains(pt)) + { + // Have to calculate directly as outside the octree + volType[pointi] = tree().shapes().getVolumeType(tree(), pt); + } + else + { + // - use cached volume type per each tree node + volType[pointi] = tree().getVolumeType(pt); + } + } + + indexedOctree::perturbTol() = oldTol; +} + + void Foam::triSurfaceMesh::setField(const labelList& values) { autoPtr fldPtr @@ -846,37 +877,6 @@ void Foam::triSurfaceMesh::getField } -void Foam::triSurfaceMesh::getVolumeType -( - const pointField& points, - List& volType -) const -{ - volType.setSize(points.size()); - - scalar oldTol = indexedOctree::perturbTol(); - indexedOctree::perturbTol() = tolerance(); - - forAll(points, pointi) - { - const point& pt = points[pointi]; - - if (!tree().bb().contains(pt)) - { - // Have to calculate directly as outside the octree - volType[pointi] = tree().shapes().getVolumeType(tree(), pt); - } - else - { - // - use cached volume type per each tree node - volType[pointi] = tree().getVolumeType(pt); - } - } - - indexedOctree::perturbTol() = oldTol; -} - - bool Foam::triSurfaceMesh::writeObject ( IOstream::streamFormat fmt, diff --git a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H index 338da94ac..2831a8cff 100644 --- a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H +++ b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -51,6 +51,8 @@ SourceFiles #include "EdgeMap.H" #include "triSurface.H" #include "triSurfaceRegionSearch.H" +#include "triSurfaceFieldsFwd.H" +#include "pointIndexHitList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -131,6 +133,30 @@ class triSurfaceMesh DynamicList& hits ); + void drawHitProblem + ( + const label fi, + const point& start, + const point& p, + const point& end, + const pointIndexHitList& hitInfo + ) const; + + void processHit + ( + scalar& internalCloseness, + scalar& externalCloseness, + const scalar internalToleranceCosAngle, + const scalar externalToleranceCosAngle, + const label fi, + const point& start, + const point& p, + const point& end, + const vector& normal, + const vectorField& normals, + const pointIndexHitList& hitInfo + ) const; + //- Disallow default bitwise copy construct triSurfaceMesh(const triSurfaceMesh&); @@ -290,6 +316,18 @@ public: // indices) get the specified field. Misses do not get set. virtual void getField(const List&, labelList&) const; + Pair> extractCloseness + ( + const scalar internalAngleTolerance = 80, + const scalar externalAngleTolerance = 10 + ) const; + + Pair> extractPointCloseness + ( + const scalar internalAngleTolerance = 80, + const scalar externalAngleTolerance = 10 + ) const; + // regIOobject implementation diff --git a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C index a1d525210..ed4d8edf9 100644 --- a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C +++ b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C @@ -1499,4 +1499,167 @@ void Foam::surfaceFeatures::operator=(const surfaceFeatures& rhs) } +// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * // + +Foam::surfaceFeatures::edgeStatus Foam::checkNonManifoldEdge +( + const triSurface& surf, + const scalar tol, + const scalar includedAngle, + const label edgei +) +{ + const edge& e = surf.edges()[edgei]; + const labelList& eFaces = surf.edgeFaces()[edgei]; + + // Bin according to normal + + DynamicList normals(2); + DynamicList bins(2); + + forAll(eFaces, eFacei) + { + const Foam::vector& n = surf.faceNormals()[eFaces[eFacei]]; + + // Find the normal in normals + label index = -1; + forAll(normals, normalI) + { + if (mag(n&normals[normalI]) > (1-tol)) + { + index = normalI; + break; + } + } + + if (index != -1) + { + bins[index].append(eFacei); + } + else if (normals.size() >= 2) + { + // Would be third normal. Mark as feature. + //Pout<< "** at edge:" << surf.localPoints()[e[0]] + // << surf.localPoints()[e[1]] + // << " have normals:" << normals + // << " and " << n << endl; + return surfaceFeatures::REGION; + } + else + { + normals.append(n); + bins.append(labelList(1, eFacei)); + } + } + + // Check resulting number of bins + if (bins.size() == 1) + { + // Note: should check here whether they are two sets of faces + // that are planar or indeed 4 faces al coming together at an edge. + //Pout<< "** at edge:" + // << surf.localPoints()[e[0]] + // << surf.localPoints()[e[1]] + // << " have single normal:" << normals[0] + // << endl; + return surfaceFeatures::NONE; + } + else + { + // Two bins. Check if normals make an angle + + //Pout<< "** at edge:" + // << surf.localPoints()[e[0]] + // << surf.localPoints()[e[1]] << nl + // << " normals:" << normals << nl + // << " bins :" << bins << nl + // << endl; + + if (includedAngle >= 0) + { + scalar minCos = Foam::cos(degToRad(180.0 - includedAngle)); + + forAll(eFaces, i) + { + const Foam::vector& ni = surf.faceNormals()[eFaces[i]]; + for (label j=i+1; j 0) + { + regionAndNormal[i] = t.region()+1; + } + else if (dir == 0) + { + FatalErrorInFunction + << exit(FatalError); + } + else + { + regionAndNormal[i] = -(t.region()+1); + } + } + + // 2. check against bin1 + const labelList& bin1 = bins[1]; + labelList regionAndNormal1(bin1.size()); + forAll(bin1, i) + { + const labelledTri& t = surf.localFaces()[eFaces[bin1[i]]]; + int dir = t.edgeDirection(e); + + label myRegionAndNormal; + if (dir > 0) + { + myRegionAndNormal = t.region()+1; + } + else + { + myRegionAndNormal = -(t.region()+1); + } + + regionAndNormal1[i] = myRegionAndNormal; + + label index = findIndex(regionAndNormal, -myRegionAndNormal); + if (index == -1) + { + // Not found. + //Pout<< "cannot find region " << myRegionAndNormal + // << " in regions " << regionAndNormal << endl; + + return surfaceFeatures::REGION; + } + } + + return surfaceFeatures::NONE; + } +} + + // ************************************************************************* // diff --git a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H index 3eb4e6f02..1cb73687b 100644 --- a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H +++ b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -32,13 +32,11 @@ Description externalStart_ .. internalStart_-1 : external edges internalStart_ .. size-1 : internal edges - NOTE: angle is included angle, not feature angle and is in degrees. The included angle is the smallest angle between two planes. For coplanar faces it is 180, for straight angles it is 90. To pick up straight edges only use included angle of 91 degrees - SourceFiles surfaceFeatures.C @@ -418,11 +416,22 @@ public: // Member Operators void operator=(const surfaceFeatures&); - - }; +//- Divide into multiple normal bins +// - return REGION if != 2 normals +// - return REGION if 2 normals that make feature angle +// - otherwise return NONE and set normals,bins +surfaceFeatures::edgeStatus checkNonManifoldEdge +( + const triSurface& surf, + const scalar tol, + const scalar includedAngle, + const label edgei +); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam