zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: cvMesh: Add vectorTools functions, a test folder and update cvMesh

Browse Source
This commit is contained in:
laurence
2012-01-12 16:03:48 +00:00
parent 6f3027e221
commit 674abd8ecd
7 changed files with 350 additions and 88 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
applications/test/vectorTools/Make/files Normal file
View File

@ -0,0 +1,3 @@
Test-vectorTools.C
EXE = $(FOAM_USER_APPBIN)/Test-vectorTools

1
applications/test/vectorTools/Make/options Normal file
View File

@ -0,0 +1 @@
EXE_INC = -I$(FOAM_APP)/utilities/mesh/generation/cvMesh/vectorTools

71
applications/test/vectorTools/Test-vectorTools.C Normal file
View File

@ -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;
}

58
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
View File

@ -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<pointIndexHit>(1, info),
List<pointIndexHit>(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<pointIndexHit>(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,

155
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C
View File

@ -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<<endl;
// Todo,needed later but want to get rid of this.
const Foam::point concaveEdgeLocalFeatPt = featPt + ppDist*concaveEdgeDir;
@ -246,6 +238,9 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
forAll(concaveEdgeNormals, edgeNormalI)
{
bool convexEdgeA = false;
bool convexEdgeB = false;
forAll(convexEdgeANormals, edgeAnormalI)
{
const vector& concaveNormal
@ -253,21 +248,9 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
const vector& convexNormal
= normals[convexEdgeANormals[edgeAnormalI]];
const scalar orientation = concaveNormal & convexNormal;
if (orientation <= 0.0)
if (areParallel(concaveNormal, convexNormal))
{
convexEdgePlaneCNormal = convexNormal;
plane planeC(featPt, convexEdgePlaneCNormal);
externalPtD =
internalPtA
+ 2.0*planeC.distance(internalPtA)*convexEdgePlaneCNormal;
pts.append(externalPtD);
indices.append(0);
types.append(-2);
convexEdgeA = true;
}
}
@ -278,21 +261,84 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
const vector& convexNormal
= normals[convexEdgeBNormals[edgeBnormalI]];
const scalar orientation = concaveNormal & convexNormal;
if (orientation <= 0.0)
// Need a looser tolerance, because sometimes adjacent triangles
// on the same surface will be slightly out of alignment.
if (areParallel(concaveNormal, convexNormal))
{
convexEdgePlaneDNormal = convexNormal;
convexEdgeB = true;
}
}
plane planeD(featPt, convexEdgePlaneDNormal);
if (convexEdgeA && convexEdgeB)
{
FatalErrorIn
(
"Foam::conformalVoronoiMesh"
"::createSpecialisedFeaturePoint"
)
<< "Both convex edges share the concave edges normal!"
<< exit(FatalError);
}
externalPtE =
internalPtB
+ 2.0*planeD.distance(internalPtB)*convexEdgePlaneDNormal;
if (convexEdgeA)
{
forAll(convexEdgeANormals, edgeAnormalI)
{
const vector& concaveNormal
= normals[concaveEdgeNormals[edgeNormalI]];
const vector& convexNormal
= normals[convexEdgeANormals[edgeAnormalI]];
pts.append(externalPtE);
indices.append(0);
types.append(-2);
if
(
areObtuse(concaveNormal, convexNormal)
|| areOrthogonal(concaveNormal, convexNormal)
)
{
convexEdgePlaneCNormal = convexNormal;
plane planeC(featPt, convexEdgePlaneCNormal);
externalPtD =
internalPtA
+ 2.0*planeC.distance(internalPtA)
*convexEdgePlaneCNormal;
pts.append(externalPtD);
indices.append(0);
types.append(-2);
}
}
}
if (convexEdgeB)
{
forAll(convexEdgeBNormals, edgeBnormalI)
{
const vector& concaveNormal
= normals[concaveEdgeNormals[edgeNormalI]];
const vector& convexNormal
= normals[convexEdgeBNormals[edgeBnormalI]];
if
(
areObtuse(concaveNormal, convexNormal)
|| areOrthogonal(concaveNormal, convexNormal)
)
{
convexEdgePlaneDNormal = convexNormal;
plane planeD(featPt, convexEdgePlaneDNormal);
externalPtE =
internalPtB
+ 2.0*planeD.distance(internalPtB)
*convexEdgePlaneDNormal;
pts.append(externalPtE);
indices.append(0);
types.append(-2);
}
}
}
}
@ -303,19 +349,10 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
const scalar totalAngle = radToDeg
(
constant::mathematical::pi +
acos(mag(concaveEdgePlaneANormal & concaveEdgePlaneBNormal))
constant::mathematical::pi
+ radAngleBetween(concaveEdgePlaneANormal, concaveEdgePlaneBNormal)
);
if (debug)
{
Info<< internalPtA << " " << internalPtB << nl
<< externalPtD << " " << externalPtE << endl;
Info<< "normals: " << nl
<< normals[concaveEdgeNormals[0]] << nl
<< normals[concaveEdgeNormals[1]] << endl;
}
if (totalAngle > 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]);
}
}
}

9
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C
View File

@ -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;

141
applications/utilities/mesh/generation/cvMesh/vectorTools/vectorTools.H Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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 <typename T>
bool areParallel
(
const Vector<T>& a,
const Vector<T>& 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 <typename T>
bool areOrthogonal
(
const Vector<T>& a,
const Vector<T>& 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 <typename T>
bool areAcute
(
const Vector<T>& a,
const Vector<T>& b
)
{
return ((a & b) > 0) ? true : false;
}
//- Test if angle between a and b is obtuse: a.b < 0
template <typename T>
bool areObtuse
(
const Vector<T>& a,
const Vector<T>& b
)
{
return ((a & b) < 0) ? true : false;
}
//- Calculate angle between a and b in radians
template <typename T>
T radAngleBetween
(
const Vector<T>& a,
const Vector<T>& b,
const T& tolerance = SMALL
)
{
return Foam::acos( (a & b)/(mag(a)*mag(b) + tolerance) );
}
//- Calculate angle between a and b in degrees
template <typename T>
T degAngleBetween
(
const Vector<T>& a,
const Vector<T>& b,
const T& tolerance = SMALL
)
{
return Foam::radToDeg(radAngleBetween(a, b, tolerance));
}
} // End namespace vectorTools
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //