From c618e6a9d3940521f4dd81adb104c5c2414d3dbc Mon Sep 17 00:00:00 2001 From: laurence Date: Fri, 13 Jan 2012 12:11:22 +0000 Subject: [PATCH] BUG: cvMesh: added tolerance for deciding whether vectors are parallel. There was a bug in feature point handling. normals that were effectively parallel were not being picked up as being parallel, so a tolerance has been added as a static const. --- .../cvMesh/conformalVoronoiMesh/Make/options | 3 +- .../conformalVoronoiMesh.C | 2 + .../conformalVoronoiMesh.H | 4 +- .../conformalVoronoiMeshConformToSurface.C | 89 ++++++++++++++----- ...alVoronoiMeshFeaturePointSpecialisations.C | 22 +++-- .../cvMesh/vectorTools/vectorTools.H | 11 +-- 6 files changed, 87 insertions(+), 44 deletions(-) diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options index 711e4ae773..1a1902b559 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options @@ -14,7 +14,8 @@ EXE_INC = \ -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \ -I$(LIB_SRC)/edgeMesh/lnInclude \ -I$(LIB_SRC)/dynamicMesh/lnInclude \ - -I$(LIB_SRC)/triSurface/lnInclude + -I$(LIB_SRC)/triSurface/lnInclude \ + -I../vectorTools EXE_LIBS = \ -lmeshTools \ diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C index c93d3236c8..dc19b47b69 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C @@ -39,6 +39,8 @@ defineTypeNameAndDebug(conformalVoronoiMesh, 0); } +const Foam::scalar Foam::conformalVoronoiMesh::tolParallel = 1e-3; + // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H index 7653bdc3d9..1946b56d7b 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H @@ -121,6 +121,8 @@ public: private: + static const scalar tolParallel; + // Private data //- The time registry of the application @@ -658,7 +660,7 @@ private: //- edge conformation location bool nearFeatureEdgeLocation ( - const Foam::point& pt, + pointIndexHit& pHit, DynamicList& newEdgeLocations, pointField& existingEdgeLocations, autoPtr >& edgeLocationTree diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C index 9fad318a8b..8dcfff585a 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C @@ -1821,40 +1821,45 @@ void Foam::conformalVoronoiMesh::limitDisplacement bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation ( - const Foam::point& pt, + pointIndexHit& pHit, DynamicList& newEdgeLocations, pointField& existingEdgeLocations, autoPtr >& edgeLocationTree ) const { + const Foam::point pt = pHit.hitPoint(); + scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt); // 0.01 and 1000 determined from speed tests, varying the indexedOctree // rebuild frequency and balance of additions to queries. - if - ( - newEdgeLocations.size() - >= max(0.01*existingEdgeLocations.size(), 1000) - ) - { +// if +// ( +// newEdgeLocations.size() +// >= max(0.01*existingEdgeLocations.size(), 1000) +// ) +// { + const label originalSize = existingEdgeLocations.size(); + existingEdgeLocations.append(newEdgeLocations); buildEdgeLocationTree(edgeLocationTree, existingEdgeLocations); newEdgeLocations.clear(); - } - else - { - // Search for the nearest point in newEdgeLocations. - // Searching here first, because the intention is that the value of - // newEdgeLocationsSizeLimit should make this faster by design. - - if (min(magSqr(newEdgeLocations - pt)) <= exclusionRangeSqr) - { - return true; - } - } +// } +// else +// { +// // Search for the nearest point in newEdgeLocations. +// // Searching here first, because the intention is that the value of +// // newEdgeLocationsSizeLimit should make this faster by design. +// +// if (min(magSqr(newEdgeLocations - pt)) <= exclusionRangeSqr) +// { +// // Average the points... +// return true; +// } +// } // Searching for the nearest point in existingEdgeLocations using the // indexedOctree @@ -1862,6 +1867,44 @@ bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation pointIndexHit info = edgeLocationTree().findNearest(pt, exclusionRangeSqr); // Average the points... +// if (info.hit()) +// { +// Foam::point newPt = 0.5*(info.hitPoint() + pt); +// +// pHit.setPoint(newPt); +// +// boolList toRemove(existingEdgeLocations.size(), false); +// +// forAll(existingEdgeLocations, pI) +// { +// const Foam::point& existingPoint = existingEdgeLocations[pI]; +// +// if (existingPoint == info.hitPoint()) +// { +// toRemove[pI] = true; +// //Info<< "Replace " << pt << " and " << info.hitPoint() << " with " << newPt << endl; +// } +// } +// +// pointField newExistingEdgeLocations(existingEdgeLocations.size()); +// +// label count = 0; +// forAll(existingEdgeLocations, pI) +// { +// if (toRemove[pI] == false) +// { +// newExistingEdgeLocations[count++] = existingEdgeLocations[pI]; +// } +// } +// +// newExistingEdgeLocations.resize(count); +// +// existingEdgeLocations = newExistingEdgeLocations; +// +// existingEdgeLocations.append(newPt); +// +// return !info.hit(); +// } return info.hit(); } @@ -2048,7 +2091,7 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits ( !nearFeatureEdgeLocation ( - edHit.hitPoint(), + edHit, newEdgeLocations, existingEdgeLocations, edgeLocationTree @@ -2114,9 +2157,9 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits continue; } - const pointIndexHit& edHit = edHits[i]; + pointIndexHit& edHit = edHits[i]; - const label featureHit = featuresHit[i]; + label featureHit = featuresHit[i]; if (edHit.hit()) { @@ -2144,7 +2187,7 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits ( !nearFeatureEdgeLocation ( - edHit.hitPoint(), + edHit, newEdgeLocations, existingEdgeLocations, edgeLocationTree diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C index 688113916f..28237dc4a8 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C @@ -248,7 +248,7 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint const vector& convexNormal = normals[convexEdgeANormals[edgeAnormalI]]; - if (areParallel(concaveNormal, convexNormal)) + if (areParallel(concaveNormal, convexNormal, tolParallel)) { convexEdgeA = true; } @@ -263,7 +263,7 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint // Need a looser tolerance, because sometimes adjacent triangles // on the same surface will be slightly out of alignment. - if (areParallel(concaveNormal, convexNormal)) + if (areParallel(concaveNormal, convexNormal, tolParallel)) { convexEdgeB = true; } @@ -291,8 +291,7 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint if ( - areObtuse(concaveNormal, convexNormal) - || areOrthogonal(concaveNormal, convexNormal) + !areParallel(concaveNormal, convexNormal, tolParallel) ) { convexEdgePlaneCNormal = convexNormal; @@ -322,8 +321,7 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint if ( - areObtuse(concaveNormal, convexNormal) - || areOrthogonal(concaveNormal, convexNormal) + !areParallel(concaveNormal, convexNormal, tolParallel) ) { convexEdgePlaneDNormal = convexNormal; @@ -400,14 +398,14 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint pts.append(externalPtG); indices.append(0); types.append(-1); + } - if (debug) + if (debug) + { + for (label ptI = initialNumOfPoints; ptI < pts.size(); ++ptI) { - for (label ptI = initialNumOfPoints; ptI < pts.size(); ++ptI) - { - Info<< "Point " << ptI << " : "; - meshTools::writeOBJ(Info, pts[ptI]); - } + Info<< "Point " << ptI << " : "; + meshTools::writeOBJ(Info, pts[ptI]); } } diff --git a/applications/utilities/mesh/generation/cvMesh/vectorTools/vectorTools.H b/applications/utilities/mesh/generation/cvMesh/vectorTools/vectorTools.H index 1228074ffe..91ffae198d 100644 --- a/applications/utilities/mesh/generation/cvMesh/vectorTools/vectorTools.H +++ b/applications/utilities/mesh/generation/cvMesh/vectorTools/vectorTools.H @@ -21,11 +21,6 @@ License You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see . - As a special exception, you have permission to link this program with the - CGAL library and distribute executables, as long as you follow the - requirements of the GNU GPL in regard to all of the software in the - executable aside from CGAL. - Class Foam::vectorTools @@ -54,8 +49,9 @@ namespace Foam //- Collection of functions for testing relationships between two vectors. namespace vectorTools { - //- Test if a and b are parallel: a.b = 1 + // Uses the cross product, so the tolerance is proportional to + // the sine of the angle between a and b in radians template bool areParallel ( @@ -71,6 +67,8 @@ namespace vectorTools } //- Test if a and b are orthogonal: a.b = 0 + // Uses the dot product, so the tolerance is proportional to + // the cosine of the angle between a and b in radians template bool areOrthogonal ( @@ -127,7 +125,6 @@ namespace vectorTools { return Foam::radToDeg(radAngleBetween(a, b, tolerance)); } - } // End namespace vectorTools // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //