Commit e3f20df1 authored by mattijs's avatar mattijs
Browse files

ENH: checkMesh : wedge & empty checking improved

- aspect ratio and cellDeterminant do not use 3rd direction
- wedges are properly check for having opposite one
parent c7bc9216
#include "checkGeometry.H"
#include "polyMesh.H"
#include "globalMeshData.H"
#include "cellSet.H"
#include "faceSet.H"
#include "pointSet.H"
#include "EdgeMap.H"
#include "wedgePolyPatch.H"
#include "unitConversion.H"
// Find wedge with opposite orientation. Note: does not actually check that
// it is opposite, only that it has opposite normal and same axis
Foam::label Foam::findOppositeWedge
(
const polyMesh& mesh,
const wedgePolyPatch& wpp
)
{
const polyBoundaryMesh& patches = mesh.boundaryMesh();
scalar wppCosAngle = wpp.centreNormal()&wpp.patchNormal();
forAll(patches, patchI)
{
if
(
patchI != wpp.index()
&& patches[patchI].size()
&& isA<wedgePolyPatch>(patches[patchI])
)
{
const wedgePolyPatch& pp = refCast<const wedgePolyPatch>
(
patches[patchI]
);
// Calculate (cos of) angle to wpp (not pp!) centre normal
scalar ppCosAngle = wpp.centreNormal()&pp.patchNormal();
if
(
pp.size() == wpp.size()
&& mag(pp.axis() & wpp.axis()) >= (1-1E-3)
&& mag(ppCosAngle - wppCosAngle) >= 1E-3
)
{
return patchI;
}
}
}
return -1;
}
bool Foam::checkWedges
(
const polyMesh& mesh,
const bool report,
const Vector<label>& directions,
labelHashSet* setPtr
)
{
// To mark edges without calculating edge addressing
EdgeMap<label> edgesInError;
const pointField& p = mesh.points();
const faceList& fcs = mesh.faces();
const polyBoundaryMesh& patches = mesh.boundaryMesh();
forAll(patches, patchI)
{
if (patches[patchI].size() && isA<wedgePolyPatch>(patches[patchI]))
{
const wedgePolyPatch& pp = refCast<const wedgePolyPatch>
(
patches[patchI]
);
scalar wedgeAngle = acos(pp.centreNormal()&pp.patchNormal());
if (report)
{
Info<< " Wedge " << pp.name() << " with angle "
<< radToDeg(wedgeAngle) << " degrees"
<< endl;
}
// Find opposite
label oppositePatchI = findOppositeWedge(mesh, pp);
if (oppositePatchI == -1)
{
if (report)
{
Info<< " ***Cannot find opposite wedge for wedge "
<< pp.name() << endl;
}
return true;
}
const wedgePolyPatch& opp = refCast<const wedgePolyPatch>
(
patches[oppositePatchI]
);
if (mag(opp.axis() & pp.axis()) < (1-1E-3))
{
if (report)
{
Info<< " ***Wedges do not have the same axis."
<< " Encountered " << pp.axis()
<< " on patch " << pp.name()
<< " which differs from " << opp.axis()
<< " on opposite wedge patch" << opp.axis()
<< endl;
}
return true;
}
// Mark edges on wedgePatches
forAll(pp, i)
{
const face& f = pp[i];
forAll(f, fp)
{
label p0 = f[fp];
label p1 = f.nextLabel(fp);
edgesInError.insert(edge(p0, p1), -1); // non-error value
}
}
// Check that wedge patch is flat
const point& p0 = p[pp.meshPoints()[0]];
forAll(pp.meshPoints(), i)
{
const point& pt = p[pp.meshPoints()[i]];
scalar d = mag((pt-p0) & pp.patchNormal());
if (d > sqrt(SMALL))
{
if (report)
{
Info<< " ***Wedge patch " << pp.name() << " not planar."
<< " Point " << pt << " is not in patch plane by "
<< d << " meter."
<< endl;
}
return true;
}
}
}
}
// Check all non-wedge faces
label nEdgesInError = 0;
forAll(fcs, faceI)
{
const face& f = fcs[faceI];
forAll(f, fp)
{
label p0 = f[fp];
label p1 = f.nextLabel(fp);
if (p0 < p1)
{
vector d(p[p1]-p[p0]);
scalar magD = mag(d);
if (magD > ROOTVSMALL)
{
d /= magD;
// Check how many empty directions are used by the edge.
label nEmptyDirs = 0;
label nNonEmptyDirs = 0;
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
if (mag(d[cmpt]) > 1e-6)
{
if (directions[cmpt] == 0)
{
nEmptyDirs++;
}
else
{
nNonEmptyDirs++;
}
}
}
if (nEmptyDirs == 0)
{
// Purely in ok directions.
}
else if (nEmptyDirs == 1)
{
// Ok if purely in empty directions.
if (nNonEmptyDirs > 0)
{
if (edgesInError.insert(edge(p0, p1), faceI))
{
nEdgesInError++;
}
}
}
else if (nEmptyDirs > 1)
{
// Always an error
if (edgesInError.insert(edge(p0, p1), faceI))
{
nEdgesInError++;
}
}
}
}
}
}
label nErrorEdges = returnReduce(nEdgesInError, sumOp<label>());
if (nErrorEdges > 0)
{
if (report)
{
Info<< " ***Number of edges not aligned with or perpendicular to "
<< "non-empty directions: " << nErrorEdges << endl;
}
if (setPtr)
{
setPtr->resize(2*nEdgesInError);
forAllConstIter(EdgeMap<label>, edgesInError, iter)
{
if (iter() >= 0)
{
setPtr->insert(iter.key()[0]);
setPtr->insert(iter.key()[1]);
}
}
}
return true;
}
else
{
if (report)
{
Info<< " All edges aligned with or perpendicular to "
<< "non-empty directions." << endl;
}
return false;
}
}
Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
{
......@@ -33,7 +288,17 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
{
pointSet nonAlignedPoints(mesh, "nonAlignedEdges", mesh.nPoints()/100);
if (mesh.checkEdgeAlignment(true, validDirs, &nonAlignedPoints))
if
(
(
validDirs != solDirs
&& checkWedges(mesh, true, validDirs, &nonAlignedPoints)
)
|| (
validDirs == solDirs
&& mesh.checkEdgeAlignment(true, validDirs, &nonAlignedPoints)
)
)
{
noFailedChecks++;
label nNonAligned = returnReduce
......@@ -58,7 +323,16 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
{
cellSet cells(mesh, "nonClosedCells", mesh.nCells()/100+1);
cellSet aspectCells(mesh, "highAspectRatioCells", mesh.nCells()/100+1);
if (mesh.checkClosedCells(true, &cells, &aspectCells))
if
(
mesh.checkClosedCells
(
true,
&cells,
&aspectCells,
mesh.geometricD()
)
)
{
noFailedChecks++;
......@@ -278,7 +552,7 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
if (allGeometry)
{
cellSet cells(mesh, "underdeterminedCells", mesh.nCells()/100);
if (mesh.checkCellDeterminant(true, &cells))
if (mesh.checkCellDeterminant(true, &cells, mesh.geometricD()))
{
noFailedChecks++;
......
#include "label.H"
#include "HashSet.H"
#include "labelVector.H"
namespace Foam
{
class polyMesh;
class wedgePolyPatch;
label findOppositeWedge(const polyMesh&, const wedgePolyPatch&);
bool checkWedges
(
const polyMesh&,
const bool report,
const Vector<label>&,
labelHashSet*
);
label checkGeometry(const polyMesh& mesh, const bool allGeometry);
}
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -544,7 +544,8 @@ public:
(
const bool report = false,
labelHashSet* setPtr = NULL,
labelHashSet* highAspectSetPtr = NULL
labelHashSet* highAspectSetPtr = NULL,
const Vector<label>& solutionD = Vector<label>::one
) const;
//- Check for negative face areas
......@@ -645,7 +646,8 @@ public:
bool checkCellDeterminant
(
const bool report = false,
labelHashSet* setPtr = NULL
labelHashSet* setPtr = NULL,
const Vector<label>& solutionD = Vector<label>::one
) const;
//- Check for concave cells by the planes of faces
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -98,13 +98,15 @@ bool Foam::primitiveMesh::checkClosedCells
(
const bool report,
labelHashSet* setPtr,
labelHashSet* aspectSetPtr
labelHashSet* aspectSetPtr,
const Vector<label>& meshD
) const
{
if (debug)
{
Info<< "bool primitiveMesh::checkClosedCells("
<< "const bool, labelHashSet*, labelHashSet*) const: "
<< "const bool, labelHashSet*, labelHashSet*"
<< ", const Vector<label>&) const: "
<< "checking whether cells are closed" << endl;
}
......@@ -171,6 +173,16 @@ bool Foam::primitiveMesh::checkClosedCells
const scalarField& vols = cellVolumes();
label nDims = 0;
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (meshD[dir] == 1)
{
nDims++;
}
}
// Check the sums
forAll(sumClosed, cellI)
{
......@@ -200,12 +212,26 @@ bool Foam::primitiveMesh::checkClosedCells
// Calculate the aspect ration as the maximum of Cartesian component
// aspect ratio to the total area hydraulic area aspect ratio
scalar aspectRatio = max
(
cmptMax(sumMagClosed[cellI])
/(cmptMin(sumMagClosed[cellI]) + VSMALL),
1.0/6.0*cmptSum(sumMagClosed[cellI])/pow(vols[cellI], 2.0/3.0)
);
scalar minCmpt = VGREAT;
scalar maxCmpt = -VGREAT;
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (meshD[dir] == 1)
{
minCmpt = min(minCmpt, sumMagClosed[cellI][dir]);
maxCmpt = max(maxCmpt, sumMagClosed[cellI][dir]);
}
}
scalar aspectRatio = maxCmpt/(minCmpt + VSMALL);
if (nDims == 3)
{
aspectRatio = max
(
aspectRatio,
1.0/6.0*cmptSum(sumMagClosed[cellI])/pow(vols[cellI], 2.0/3.0)
);
}
maxAspectRatio = max(maxAspectRatio, aspectRatio);
......@@ -1977,7 +2003,8 @@ bool Foam::primitiveMesh::checkFaceFaces
bool Foam::primitiveMesh::checkCellDeterminant
(
const bool report, // report,
labelHashSet* setPtr // setPtr
labelHashSet* setPtr, // setPtr
const Vector<label>& meshD
) const
{
if (debug)
......@@ -1987,6 +2014,22 @@ bool Foam::primitiveMesh::checkCellDeterminant
<< "checking for under-determined cells" << endl;
}
// Determine number of dimensions and (for 2D) missing dimension
label nDims = 0;
label twoD = -1;
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (meshD[dir] == 1)
{
nDims++;
}
else
{
twoD = dir;
}
}
const cellList& c = cells();
label nErrorCells = 0;
......@@ -1995,62 +2038,89 @@ bool Foam::primitiveMesh::checkCellDeterminant
scalar sumDet = 0;
label nSummed = 0;
forAll(c, cellI)
if (nDims == 1)
{
const labelList& curFaces = c[cellI];
minDet = 1;
sumDet = c.size()*minDet;
nSummed = c.size();
}
else
{
forAll (c, cellI)
{
const labelList& curFaces = c[cellI];
// Calculate local normalization factor
scalar avgArea = 0;
// Calculate local normalization factor
scalar avgArea = 0;
label nInternalFaces = 0;
label nInternalFaces = 0;
forAll(curFaces, i)
{
if (isInternalFace(curFaces[i]))
forAll(curFaces, i)
{
avgArea += mag(faceAreas()[curFaces[i]]);
if (isInternalFace(curFaces[i]))
{
avgArea += mag(faceAreas()[curFaces[i]]);
nInternalFaces++;
nInternalFaces++;
}
}
}
if (nInternalFaces == 0)
{
if (setPtr)
if (nInternalFaces == 0)
{
setPtr->insert(cellI);
if (setPtr)
{
setPtr->insert(cellI);
}
nErrorCells++;
}
else
{
avgArea /= nInternalFaces;
nErrorCells++;
}
else
{
avgArea /= nInternalFaces;
symmTensor areaTensor(symmTensor::zero);
symmTensor areaTensor(symmTensor::zero);
forAll(curFaces, i)
{
if (isInternalFace(curFaces[i]))
{
areaTensor += sqr(faceAreas()[curFaces[i]]/avgArea);
}
}
forAll(curFaces, i)
{
if (isInternalFace(curFaces[i]))
if (nDims == 2)
{
areaTensor += sqr(faceAreas()[curFaces[i]]/avgArea);
// Add the missing eigenvector (such that it does not
// affect the determinant)
if (twoD == 0)
{
areaTensor.xx() = 1;
}
else if (twoD == 1)
{
areaTensor.yy() = 1;
}
else
{
areaTensor.zz() = 1;
}
}
}
scalar determinant = mag(det(areaTensor));
scalar determinant = mag(det(areaTensor));
minDet = min(determinant, minDet);
sumDet += determinant;
nSummed++;
minDet = min(determinant, minDet);
sumDet += determinant;
nSummed++;