diff --git a/applications/test/vectorTools/Make/files b/applications/test/vectorTools/Make/files new file mode 100644 index 0000000000..0b30b98f8f --- /dev/null +++ b/applications/test/vectorTools/Make/files @@ -0,0 +1,3 @@ +Test-vectorTools.C + +EXE = $(FOAM_USER_APPBIN)/Test-vectorTools diff --git a/applications/test/vectorTools/Make/options b/applications/test/vectorTools/Make/options new file mode 100644 index 0000000000..9e015e6078 --- /dev/null +++ b/applications/test/vectorTools/Make/options @@ -0,0 +1 @@ +EXE_INC = -I$(FOAM_APP)/utilities/mesh/generation/cvMesh/vectorTools diff --git a/applications/test/vectorTools/Test-vectorTools.C b/applications/test/vectorTools/Test-vectorTools.C new file mode 100644 index 0000000000..0c85ac935b --- /dev/null +++ b/applications/test/vectorTools/Test-vectorTools.C @@ -0,0 +1,71 @@ +#include "vector.H" +#include "IOstreams.H" +#include "vectorTools.H" +#include "unitConversion.H" + +using namespace Foam; + + +void test(const vector& a, const vector& b, const scalar tolerance) +{ + Info<< "Vectors " << a << " and " << b + << " are (to tolerance of " << tolerance << "): "; + + if (vectorTools::areParallel(a, b, tolerance)) + Info<< " parallel "; + + if (vectorTools::areOrthogonal(a, b, tolerance)) + Info<< " orthogonal "; + + if (vectorTools::areAcute(a, b)) + Info<< " acute "; + + if (vectorTools::areObtuse(a, b)) + Info<< " obtuse "; + + Info<< ", angle = " << vectorTools::degAngleBetween(a, b); + + Info<< endl; +} + + +int main() +{ + vector a(1.0, 1.0, 1.0); + vector b(2.0, 2.0, 2.0); + + test(a, b, 0.0); + test(a, b, VSMALL); + test(a, b, SMALL); + test(a, b, 1e-3); + test(a, b, 1e-1); + + a = vector(1,0,0); + b = vector(0,2,0); + + test(a, b, 0.0); + test(a, b, VSMALL); + test(a, b, SMALL); + test(a, b, 1e-3); + test(a, b, 1e-1); + + a = vector(1,0,0); + b = vector(-1,0,0); + + test(a, b, 0.0); + test(a, b, VSMALL); + test(a, b, SMALL); + test(a, b, 1e-3); + test(a, b, 1e-1); + + a = vector(1,0,0); + b = vector(-1,2,0); + + test(a, b, 0.0); + test(a, b, VSMALL); + test(a, b, SMALL); + test(a, b, 1e-3); + test(a, b, 1e-1); + + return 0; +} diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C index ba98a06ac7..9fad318a8b 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C @@ -25,6 +25,9 @@ License #include "conformalVoronoiMesh.H" #include "backgroundMeshDecomposition.H" +#include "vectorTools.H" + +using namespace Foam::vectorTools; // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // @@ -644,31 +647,40 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections continue; } - const Foam::point dE0 = topoint(dual(fit->first)); + Foam::point dE0 = topoint(dual(fit->first)); - const Foam::point dE1 - = topoint(dual(fit->first->neighbor(fit->second))); + Foam::point dE1 = topoint(dual(fit->first->neighbor(fit->second))); - pointIndexHit info; - label hitSurface = -1; + pointIndexHit infoIntersection; + label hitSurfaceIntersection = -1; - geometryToConformTo_.findSurfaceAnyIntersection + if (Pstream::parRun()) + { + bool inProc = clipLineToProc(dE0, dE1); + + if (!inProc) + { + continue; + } + } + + geometryToConformTo_.findSurfaceNearestIntersection ( dE0, dE1, - info, - hitSurface + infoIntersection, + hitSurfaceIntersection ); - if (info.hit()) + if (infoIntersection.hit()) { flagIntersection = true; vectorField norm(1); - allGeometry_[hitSurface].getNormal + allGeometry_[hitSurfaceIntersection].getNormal ( - List(1, info), + List(1, infoIntersection), norm ); @@ -676,7 +688,7 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections const Foam::point vertex = topoint(vit->point()); - const plane p(info.hitPoint(), n); + const plane p(infoIntersection.hitPoint(), n); const plane::ray r(vertex, n); @@ -684,10 +696,13 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections const Foam::point newPoint = vertex + d*n; - geometryToConformTo_.findSurfaceAnyIntersection + pointIndexHit info; + label hitSurface = -1; + + geometryToConformTo_.findSurfaceNearestIntersection ( vertex, - vertex + 1.1*(newPoint - vertex), + newPoint + SMALL*n, info, hitSurface ); @@ -1951,17 +1966,17 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits continue; } - vectorField norm(1); + vectorField norm1(1); allGeometry_[hitSurf].getNormal ( List(1, hitInfo), - norm + norm1 ); - const vector& n = norm[0]; + const vector& n1 = norm1[0]; - const plane surfPlane(hitInfo.hitPoint(), n); + const plane surfPlane(hitInfo.hitPoint(), n1); pointsUsed.append(hI); @@ -2004,18 +2019,19 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits norm2 ); + const vector& n2 = norm2[0]; + // If the normals are nearly parallel then ignore. - if (mag(mag(norm2[0] & n) - 1.0) < SMALL) + if (areParallel(n1, n2)) { continue; } - const Foam::point newEdgePoint = surfPoint + n*d; + const Foam::point newEdgePoint = surfPoint + n1*d; pointIndexHit edHit; label featureHit = -1; - // This is a bit lazy... geometryToConformTo_.findEdgeNearest ( newEdgePoint, diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C index 108db83be9..688113916f 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C @@ -24,6 +24,9 @@ License \*---------------------------------------------------------------------------*/ #include "conformalVoronoiMesh.H" +#include "vectorTools.H" + +using namespace Foam::vectorTools; // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // @@ -117,6 +120,8 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint return false; } + const label initialNumOfPoints = pts.size(); + const scalar ppDist = pointPairDistance(featPt); const vectorField& normals = feMesh.normals(); @@ -161,13 +166,6 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint const vector& concaveEdgePlaneBNormal = normals[edgeNormals[concaveEdgeI][1]]; - if (debug) - { - Info<< featPt << nl - << concaveEdgePlaneANormal << nl - << concaveEdgePlaneBNormal << endl; - } - // Intersect planes parallel to the concave edge planes offset // by ppDist and the plane defined by featPt and the edge vector. plane planeA @@ -188,12 +186,6 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint ptI ); - Info<< "pt " << featPt << nl - << "concave edge dir " << concaveEdgeDir << nl - << "convex edge dir 1 " << feMesh.edgeDirection(convexEdgesI[0], ptI) << nl - << "convex edge dir 2 " << feMesh.edgeDirection(convexEdgesI[1], ptI) << nl< cvMeshControls().maxQuadAngle()) { // Add additional mitering points @@ -358,7 +395,6 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint const Foam::point externalPtG = internalPtF - //+ 2.0*planeM.distance(internalPtF) * convexEdgesPlaneNormal; + 2.0*planeM.distance(internalPtF)*convexEdgesPlaneNormal; pts.append(externalPtG); @@ -367,20 +403,11 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint if (debug) { - Info<< nl << "# featPt " << endl; - meshTools::writeOBJ(Info, featPt); - Info<< "# internalPtA" << endl; - meshTools::writeOBJ(Info, internalPtA); - Info<< "# internalPtB" << endl; - meshTools::writeOBJ(Info, internalPtB); - Info<< "# externalPtD" << endl; - meshTools::writeOBJ(Info, externalPtD); - Info<< "# externalPtE" << endl; - meshTools::writeOBJ(Info, externalPtE); - Info<< "# internalPtF" << endl; - meshTools::writeOBJ(Info, internalPtF); - Info<< "# externalPtG" << endl; - meshTools::writeOBJ(Info, externalPtG); + for (label ptI = initialNumOfPoints; ptI < pts.size(); ++ptI) + { + Info<< "Point " << ptI << " : "; + meshTools::writeOBJ(Info, pts[ptI]); + } } } diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C index 6914793815..aa8e6722fe 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C @@ -24,6 +24,9 @@ License \*---------------------------------------------------------------------------*/ #include "conformalVoronoiMesh.H" +#include "vectorTools.H" + +using namespace Foam::vectorTools; // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // @@ -117,7 +120,7 @@ void Foam::conformalVoronoiMesh::createExternalEdgePointGroup const vector& nA = feNormals[edNormalIs[0]]; const vector& nB = feNormals[edNormalIs[1]]; - if (mag(nA & nB) > 1 - SMALL) + if (areParallel(nA, nB)) { // The normals are nearly parallel, so this is too sharp a feature to // conform to. @@ -195,7 +198,7 @@ void Foam::conformalVoronoiMesh::createInternalEdgePointGroup const vector& nA = feNormals[edNormalIs[0]]; const vector& nB = feNormals[edNormalIs[1]]; - if (mag(nA & nB) > 1 - SMALL) + if (areParallel(nA, nB)) { // The normals are nearly parallel, so this is too sharp a feature to // conform to. @@ -234,7 +237,7 @@ void Foam::conformalVoronoiMesh::createInternalEdgePointGroup Foam::point reflectedB = reflMasterPt - 2*ppDist*nB; scalar totalAngle = - radToDeg(constant::mathematical::pi + acos(nA & nB)); + radToDeg(constant::mathematical::pi + radAngleBetween(nA, nB)); // Number of quadrants the angle should be split into int nQuads = int(totalAngle/cvMeshControls().maxQuadAngle()) + 1; diff --git a/applications/utilities/mesh/generation/cvMesh/vectorTools/vectorTools.H b/applications/utilities/mesh/generation/cvMesh/vectorTools/vectorTools.H new file mode 100644 index 0000000000..1228074ffe --- /dev/null +++ b/applications/utilities/mesh/generation/cvMesh/vectorTools/vectorTools.H @@ -0,0 +1,141 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 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 . + + 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 + +Description + Functions for analysing the relationships between vectors + +SourceFiles + +\*---------------------------------------------------------------------------*/ + +#ifndef vectorTools_H +#define vectorTools_H + +#include "vector.H" +#include "unitConversion.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Namespace vectorTools Declaration +\*---------------------------------------------------------------------------*/ + +//- Collection of functions for testing relationships between two vectors. +namespace vectorTools +{ + + //- Test if a and b are parallel: a.b = 1 + template + bool areParallel + ( + const Vector& a, + const Vector& b, + const T& tolerance = SMALL + ) + { + return (mag(a ^ b) < tolerance) ? true : false; +// return ( mag( mag(a & b)/(mag(a)*mag(b)) - 1.0 ) < tolerance ) +// ? true +// : false; + } + + //- Test if a and b are orthogonal: a.b = 0 + template + bool areOrthogonal + ( + const Vector& a, + const Vector& b, + const T& tolerance = SMALL + ) + { + return (mag(a & b) < tolerance) ? true : false; + } + + //- Test if angle between a and b is acute: a.b > 0 + template + bool areAcute + ( + const Vector& a, + const Vector& b + ) + { + return ((a & b) > 0) ? true : false; + } + + //- Test if angle between a and b is obtuse: a.b < 0 + template + bool areObtuse + ( + const Vector& a, + const Vector& b + ) + { + return ((a & b) < 0) ? true : false; + } + + //- Calculate angle between a and b in radians + template + T radAngleBetween + ( + const Vector& a, + const Vector& b, + const T& tolerance = SMALL + ) + { + return Foam::acos( (a & b)/(mag(a)*mag(b) + tolerance) ); + } + + //- Calculate angle between a and b in degrees + template + T degAngleBetween + ( + const Vector& a, + const Vector& b, + const T& tolerance = SMALL + ) + { + return Foam::radToDeg(radAngleBetween(a, b, tolerance)); + } + +} // End namespace vectorTools + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* //