diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C index 3af85508be..4945eb0b09 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C @@ -97,7 +97,6 @@ Foam::tmp Foam::polyMeshTools::faceSkewness { const labelList& own = mesh.faceOwner(); const labelList& nei = mesh.faceNeighbour(); - const faceList& fcs = mesh.faces(); const polyBoundaryMesh& pbm = mesh.boundaryMesh(); tmp tskew(new scalarField(mesh.nFaces())); @@ -154,31 +153,16 @@ Foam::tmp Foam::polyMeshTools::faceSkewness { label faceI = pp.start() + i; - vector Cpf = fCtrs[faceI] - cellCtrs[own[faceI]]; + skew[faceI] = primitiveMeshTools::boundaryFaceSkewness + ( + mesh, + p, + fCtrs, + fAreas, - vector normal = fAreas[faceI]; - normal /= mag(normal) + VSMALL; - vector d = normal*(normal & Cpf); - - - // Skewness vector - vector sv = - Cpf - - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d; - vector svHat = sv/(mag(sv) + VSMALL); - - // Normalisation distance calculated as the approximate distance - // from the face centre to the edge of the face in the direction - // of the skewness - scalar fd = 0.4*mag(d) + VSMALL; - const face& f = fcs[faceI]; - forAll(f, pi) - { - fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI]))); - } - - // Normalised skewness - skew[faceI] = mag(sv)/fd; + faceI, + cellCtrs[own[faceI]] + ); } } } diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C index e94c8edbf7..241f4eef8e 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -64,6 +64,44 @@ Foam::scalar Foam::primitiveMeshTools::faceSkewness return mag(sv)/fd; } +Foam::scalar Foam::primitiveMeshTools::boundaryFaceSkewness +( + const primitiveMesh& mesh, + const pointField& p, + const vectorField& fCtrs, + const vectorField& fAreas, + + const label faceI, + const point& ownCc +) +{ + vector Cpf = fCtrs[faceI] - ownCc; + + vector normal = fAreas[faceI]; + normal /= mag(normal) + VSMALL; + vector d = normal*(normal & Cpf); + + + // Skewness vector + vector sv = + Cpf + - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d; + vector svHat = sv/(mag(sv) + VSMALL); + + // Normalisation distance calculated as the approximate distance + // from the face centre to the edge of the face in the direction + // of the skewness + scalar fd = 0.4*mag(d) + VSMALL; + const face& f = mesh.faces()[faceI]; + forAll(f, pi) + { + fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI]))); + } + + // Normalised skewness + return mag(sv)/fd; +} + Foam::scalar Foam::primitiveMeshTools::faceOrthogonality ( @@ -119,7 +157,6 @@ Foam::tmp Foam::primitiveMeshTools::faceSkewness { const labelList& own = mesh.faceOwner(); const labelList& nei = mesh.faceNeighbour(); - const faceList& fcs = mesh.faces(); tmp tskew(new scalarField(mesh.nFaces())); scalarField& skew = tskew(); @@ -145,31 +182,15 @@ Foam::tmp Foam::primitiveMeshTools::faceSkewness for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++) { - vector Cpf = fCtrs[faceI] - cellCtrs[own[faceI]]; - - vector normal = fAreas[faceI]; - normal /= mag(normal) + VSMALL; - vector d = normal*(normal & Cpf); - - - // Skewness vector - vector sv = - Cpf - - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d; - vector svHat = sv/(mag(sv) + VSMALL); - - // Normalisation distance calculated as the approximate distance - // from the face centre to the edge of the face in the direction - // of the skewness - scalar fd = 0.4*mag(d) + VSMALL; - const face& f = fcs[faceI]; - forAll(f, pi) - { - fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI]))); - } - - // Normalised skewness - skew[faceI] = mag(sv)/fd; + skew[faceI] = boundaryFaceSkewness + ( + mesh, + p, + fCtrs, + fAreas, + faceI, + cellCtrs[own[faceI]] + ); } return tskew; diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H index 5e3b2c07fd..170d089e11 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -48,32 +48,6 @@ namespace Foam class primitiveMeshTools { - -protected: - - // Protected Member Functions - - //- Skewness of single face - static scalar faceSkewness - ( - const primitiveMesh& mesh, - const pointField& p, - const vectorField& fCtrs, - const vectorField& fAreas, - - const label faceI, - const point& ownCc, - const point& neiCc - ); - - //- Orthogonality of single face - static scalar faceOrthogonality - ( - const point& ownCc, - const point& neiCc, - const vector& s - ); - public: //- Generate non-orthogonality field (internal faces only) @@ -154,6 +128,41 @@ public: const PackedBoolList& internalOrCoupledFace ); + + // Helpers: single face check + + //- Skewness of single face + static scalar faceSkewness + ( + const primitiveMesh& mesh, + const pointField& p, + const vectorField& fCtrs, + const vectorField& fAreas, + + const label faceI, + const point& ownCc, + const point& neiCc + ); + + //- Skewness of single boundary face + static scalar boundaryFaceSkewness + ( + const primitiveMesh& mesh, + const pointField& p, + const vectorField& fCtrs, + const vectorField& fAreas, + + const label faceI, + const point& ownCc + ); + + //- Orthogonality of single face + static scalar faceOrthogonality + ( + const point& ownCc, + const point& neiCc, + const vector& s + ); };