ENH: polyMeshTools: abstract boundary skewness calc

This commit is contained in:
mattijs
2013-05-31 17:38:49 +01:00
parent 487ef5ac0b
commit 4c70254450
3 changed files with 93 additions and 79 deletions

View File

@ -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::scalarField> Foam::primitiveMeshTools::faceSkewness
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const faceList& fcs = mesh.faces();
tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
scalarField& skew = tskew();
@ -145,31 +182,15 @@ Foam::tmp<Foam::scalarField> 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;