Commit b328499d authored by graham's avatar graham
Browse files

ENH: Various improvements.

  + No fatal error on triSurfaceTools::surfaceSide, commented out WarningIn.

  + Make cellSizeControlSurfaces look for a GREAT span for the nearest surface
    point.

  + Identify and limit filtering on single internal face cells in polyMesh
    quality assessment.

  + Create cellSet of remaining protruding cells after polyMesh creation.

  + Implemented wellOutside function by generalising wellInside to
    wellInOutSide.
parent 5cd9ccf8
......@@ -97,6 +97,53 @@ void deleteBox
}
void drawHitProblem
(
label fI,
const triSurface& surf,
const pointField& start,
const pointField& faceCentres,
const pointField& end,
const List<pointIndexHit>& hitInfo
)
{
Info<< nl << "# findLineAll did not hit its own face."
<< nl << "# fI " << fI
<< nl << "# start " << start[fI]
<< nl << "# f centre " << faceCentres[fI]
<< nl << "# end " << end[fI]
<< nl << "# hitInfo " << hitInfo
<< endl;
meshTools::writeOBJ(Info, start[fI]);
meshTools::writeOBJ(Info, faceCentres[fI]);
meshTools::writeOBJ(Info, end[fI]);
Info<< "l 1 2 3" << endl;
meshTools::writeOBJ(Info, surf.points()[surf[fI][0]]);
meshTools::writeOBJ(Info, surf.points()[surf[fI][1]]);
meshTools::writeOBJ(Info, surf.points()[surf[fI][2]]);
Info<< "f 4 5 6" << endl;
forAll(hitInfo, hI)
{
label hFI = hitInfo[hI].index();
meshTools::writeOBJ(Info, surf.points()[surf[hFI][0]]);
meshTools::writeOBJ(Info, surf.points()[surf[hFI][1]]);
meshTools::writeOBJ(Info, surf.points()[surf[hFI][2]]);
Info<< "f "
<< 3*hI + 7 << " "
<< 3*hI + 8 << " "
<< 3*hI + 9
<< endl;
}
}
scalarField curvature(const triSurface& surf)
{
scalarField k(surf.points().size(), 0);
......@@ -598,23 +645,24 @@ int main(int argc, char *argv[])
if (hitInfo.size() < 1)
{
Info<< nl << "# fI " << fI
<< nl << "# start " << start[fI]
<< nl << "# f centre " << faceCentres[fI]
<< nl << "# end " << end[fI]
<< endl;
drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
// Info<< nl << "# fI " << fI
// << nl << "# start " << start[fI]
// << nl << "# f centre " << faceCentres[fI]
// << nl << "# end " << end[fI]
// << endl;
meshTools::writeOBJ(Info, start[fI]);
meshTools::writeOBJ(Info, faceCentres[fI]);
meshTools::writeOBJ(Info, end[fI]);
// meshTools::writeOBJ(Info, start[fI]);
// meshTools::writeOBJ(Info, faceCentres[fI]);
// meshTools::writeOBJ(Info, end[fI]);
Info<< "l 1 2 3" << endl;
// Info<< "l 1 2 3" << endl;
meshTools::writeOBJ(Info, surf.points()[surf[fI][0]]);
meshTools::writeOBJ(Info, surf.points()[surf[fI][1]]);
meshTools::writeOBJ(Info, surf.points()[surf[fI][2]]);
// meshTools::writeOBJ(Info, surf.points()[surf[fI][0]]);
// meshTools::writeOBJ(Info, surf.points()[surf[fI][1]]);
// meshTools::writeOBJ(Info, surf.points()[surf[fI][2]]);
Info<< "f 4 5 6" << endl;
// Info<< "f 4 5 6" << endl;
FatalErrorIn(args.executable())
<< "findLineAll did not hit its own face."
......@@ -630,40 +678,42 @@ int main(int argc, char *argv[])
}
else if (hitInfo[0].index() != fI)
{
Info<< nl << "# findLineAll did not hit its own face."
<< nl << "# fI " << fI
<< nl << "# start " << start[fI]
<< nl << "# f centre " << faceCentres[fI]
<< nl << "# end " << end[fI]
<< nl << "# hitInfo " << hitInfo
<< endl;
meshTools::writeOBJ(Info, start[fI]);
meshTools::writeOBJ(Info, faceCentres[fI]);
meshTools::writeOBJ(Info, end[fI]);
Info<< "l 1 2 3" << endl;
meshTools::writeOBJ(Info, surf.points()[surf[fI][0]]);
meshTools::writeOBJ(Info, surf.points()[surf[fI][1]]);
meshTools::writeOBJ(Info, surf.points()[surf[fI][2]]);
Info<< "f 4 5 6" << endl;
forAll(hitInfo, hI)
{
label hFI = hitInfo[hI].index();
meshTools::writeOBJ(Info, surf.points()[surf[hFI][0]]);
meshTools::writeOBJ(Info, surf.points()[surf[hFI][1]]);
meshTools::writeOBJ(Info, surf.points()[surf[hFI][2]]);
Info<< "f "
<< 3*hI + 7 << " "
<< 3*hI + 8 << " "
<< 3*hI + 9
<< endl;
}
drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
// Info<< nl << "# findLineAll did not hit its own face."
// << nl << "# fI " << fI
// << nl << "# start " << start[fI]
// << nl << "# f centre " << faceCentres[fI]
// << nl << "# end " << end[fI]
// << nl << "# hitInfo " << hitInfo
// << endl;
// meshTools::writeOBJ(Info, start[fI]);
// meshTools::writeOBJ(Info, faceCentres[fI]);
// meshTools::writeOBJ(Info, end[fI]);
// Info<< "l 1 2 3" << endl;
// meshTools::writeOBJ(Info, surf.points()[surf[fI][0]]);
// meshTools::writeOBJ(Info, surf.points()[surf[fI][1]]);
// meshTools::writeOBJ(Info, surf.points()[surf[fI][2]]);
// Info<< "f 4 5 6" << endl;
// forAll(hitInfo, hI)
// {
// label hFI = hitInfo[hI].index();
// meshTools::writeOBJ(Info, surf.points()[surf[hFI][0]]);
// meshTools::writeOBJ(Info, surf.points()[surf[hFI][1]]);
// meshTools::writeOBJ(Info, surf.points()[surf[hFI][2]]);
// Info<< "f "
// << 3*hI + 7 << " "
// << 3*hI + 8 << " "
// << 3*hI + 9
// << endl;
// }
// FatalErrorIn(args.executable())
// << "findLineAll did not hit its own face."
// << exit(FatalError);
......@@ -687,40 +737,41 @@ int main(int argc, char *argv[])
if (ownHitI < 0)
{
Info<< nl << "# findLineAll did not hit its own face."
<< nl << "# fI " << fI
<< nl << "# start " << start[fI]
<< nl << "# f centre " << faceCentres[fI]
<< nl << "# end " << end[fI]
<< nl << "# hitInfo " << hitInfo
<< endl;
meshTools::writeOBJ(Info, start[fI]);
meshTools::writeOBJ(Info, faceCentres[fI]);
meshTools::writeOBJ(Info, end[fI]);
Info<< "l 1 2 3" << endl;
meshTools::writeOBJ(Info, surf.points()[surf[fI][0]]);
meshTools::writeOBJ(Info, surf.points()[surf[fI][1]]);
meshTools::writeOBJ(Info, surf.points()[surf[fI][2]]);
Info<< "f 4 5 6" << endl;
forAll(hitInfo, hI)
{
label hFI = hitInfo[hI].index();
meshTools::writeOBJ(Info, surf.points()[surf[hFI][0]]);
meshTools::writeOBJ(Info, surf.points()[surf[hFI][1]]);
meshTools::writeOBJ(Info, surf.points()[surf[hFI][2]]);
Info<< "f "
<< 3*hI + 7 << " "
<< 3*hI + 8 << " "
<< 3*hI + 9
<< endl;
}
drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
// Info<< nl << "# findLineAll did not hit its own face."
// << nl << "# fI " << fI
// << nl << "# start " << start[fI]
// << nl << "# f centre " << faceCentres[fI]
// << nl << "# end " << end[fI]
// << nl << "# hitInfo " << hitInfo
// << endl;
// meshTools::writeOBJ(Info, start[fI]);
// meshTools::writeOBJ(Info, faceCentres[fI]);
// meshTools::writeOBJ(Info, end[fI]);
// Info<< "l 1 2 3" << endl;
// meshTools::writeOBJ(Info, surf.points()[surf[fI][0]]);
// meshTools::writeOBJ(Info, surf.points()[surf[fI][1]]);
// meshTools::writeOBJ(Info, surf.points()[surf[fI][2]]);
// Info<< "f 4 5 6" << endl;
// forAll(hitInfo, hI)
// {
// label hFI = hitInfo[hI].index();
// meshTools::writeOBJ(Info, surf.points()[surf[hFI][0]]);
// meshTools::writeOBJ(Info, surf.points()[surf[hFI][1]]);
// meshTools::writeOBJ(Info, surf.points()[surf[hFI][2]]);
// Info<< "f "
// << 3*hI + 7 << " "
// << 3*hI + 8 << " "
// << 3*hI + 9
// << endl;
// }
// FatalErrorIn(args.executable())
// << "findLineAll did not hit its own face."
......
......@@ -301,7 +301,7 @@ Foam::scalar Foam::cellSizeControlSurfaces::cellSize
cvMesh_.geometryToConformTo().findSurfaceNearest
(
pt,
cvMesh_.geometryToConformTo().spanMagSqr(),
sqr(GREAT),
surfHit,
hitSurface
);
......@@ -316,11 +316,7 @@ Foam::scalar Foam::cellSizeControlSurfaces::cellSize
"bool isSurfacePoint"
") const"
)
<< "Point " << pt << "was not within "
<< cvMesh_.geometryToConformTo().spanMag()
<< " of the surface." << nl
<< "findSurfaceNearest did not find a hit across the span "
<< "of the surfaces."
<< "Point " << pt << " did not find a nearest surface point"
<< nl << exit(FatalError) << endl;
}
else
......
......@@ -64,6 +64,7 @@ SourceFiles
#include "HashSet.H"
#include "memInfo.H"
#include "point.H"
#include "cellSet.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -797,6 +798,10 @@ public:
// target cell volume, actual cell volume and equivalent
// actual cell size (cbrt(actual cell volume)).
void writeCellSizes(const fvMesh& mesh) const;
//- Find the cellSet of the boundary cells which have points that
// protrude out of the surface beyond a tolerance.
void findRemainingProtrusionSet(const fvMesh& mesh) const;
};
......
......@@ -1531,29 +1531,76 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
wrongFaces
);
label nInvalidPolyhedra = 0;
{
// Check for cells with more than 1 but fewer than 4 faces
label nInvalidPolyhedra = 0;
const cellList& cells = pMesh.cells();
const cellList& cells = pMesh.cells();
forAll(cells, cI)
{
if (cells[cI].size() < 4 && cells[cI].size() > 0)
forAll(cells, cI)
{
// Info<< "cell " << cI << " " << cells[cI]
// << " has " << cells[cI].size() << " faces."
// << endl;
if (cells[cI].size() < 4 && cells[cI].size() > 0)
{
// Info<< "cell " << cI << " " << cells[cI]
// << " has " << cells[cI].size() << " faces."
// << endl;
nInvalidPolyhedra++;
forAll(cells[cI], cFI)
{
wrongFaces.insert(cells[cI][cFI]);
}
}
}
Info<< " cells with more than 1 but fewer than 4 faces : "
<< nInvalidPolyhedra << endl;
// Check for cells with one internal face only
labelList nInternalFaces(pMesh.nCells(), 0);
for (label fI = 0; fI < pMesh.nInternalFaces(); fI++)
{
nInternalFaces[pMesh.faceOwner()[fI]]++;
nInternalFaces[pMesh.faceNeighbour()[fI]]++;
}
nInvalidPolyhedra++;
const polyBoundaryMesh& patches = pMesh.boundaryMesh();
forAll(cells[cI], cFI)
forAll(patches, patchI)
{
if (patches[patchI].coupled())
{
const labelUList& owners = patches[patchI].faceCells();
forAll(owners, i)
{
nInternalFaces[owners[i]]++;
}
}
}
label oneInternalFaceCells = 0;
forAll(nInternalFaces, cI)
{
if (nInternalFaces[cI] <= 1)
{
wrongFaces.insert(cells[cI][cFI]);
oneInternalFaceCells++;
forAll(cells[cI], cFI)
{
wrongFaces.insert(cells[cI][cFI]);
}
}
}
Info<< " cells with with zero or one non-boundary face : "
<< oneInternalFaceCells << endl;
}
Info<< "Cells with more than 1 but fewer than 4 faces : "
<< nInvalidPolyhedra << endl;
PackedBoolList ptToBeLimited(pts.size(), false);
......
......@@ -328,6 +328,8 @@ void Foam::conformalVoronoiMesh::writeMesh
// cellCs.write();
writeCellSizes(mesh);
findRemainingProtrusionSet(mesh);
}
......@@ -512,4 +514,69 @@ void Foam::conformalVoronoiMesh::writeCellSizes
}
void Foam::conformalVoronoiMesh::findRemainingProtrusionSet
(
const fvMesh& mesh
) const
{
timeCheck("Start findRemainingProtrusionSet");
const polyBoundaryMesh& patches = mesh.boundaryMesh();
labelHashSet protrudingBoundaryPoints;
forAll(patches, patchI)
{
const labelList& patchLocalPtIs = patches[patchI].boundaryPoints();
forAll(patchLocalPtIs, ppI)
{
label meshPtI = patches[patchI].meshPoints()[patchLocalPtIs[ppI]];
const Foam::point& pt = mesh.points()[meshPtI];
if
(
geometryToConformTo_.wellOutside
(
pt,
sqr(2.0*maxSurfaceProtrusion(pt))
)
)
{
protrudingBoundaryPoints.insert(meshPtI);
}
}
}
cellSet protrudingCells
(
mesh,
"cvMesh_remainingProtrusions",
mesh.nCells()/1000
);
forAllConstIter(labelHashSet, protrudingBoundaryPoints, iter)
{
const label pointI = iter.key();
const labelList& pCells = mesh.pointCells()[pointI];
forAll(pCells, pCellI)
{
protrudingCells.insert(pCells[pCellI]);
}
}
if (!protrudingCells.empty())
{
Info<< nl << "Found " << protrudingCells.size()
<< " cells protruding from the surface, writing cellSet "
<< protrudingCells.name()
<< endl;
protrudingCells.write();
}
}
// ************************************************************************* //
......@@ -329,10 +329,11 @@ bool Foam::conformationSurfaces::outside
}
Foam::Field<bool> Foam::conformationSurfaces::wellInside
Foam::Field<bool> Foam::conformationSurfaces::wellInOutSide
(
const pointField& samplePts,
const scalarField& testDistSqr
const scalarField& testDistSqr,
bool testForInside
) const
{
List<List<searchableSurface::volumeType> > surfaceVolumeTests
......@@ -360,7 +361,7 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInside
// reference value and if the points are inside the surface by a given
// distanceSquared
Field<bool> insidePoints(samplePts.size(), true);
Field<bool> inOutSidePoint(samplePts.size(), true);
//Check if the points are inside the surface by the given distance squared
......@@ -384,23 +385,39 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInside
if (pHit.hit())
{
insidePoints[i] = false;
// If the point is within range of the surface, then it can't be
// well (in|out)side
inOutSidePoint[i] = false;
continue;
}
forAll(surfaces_, s)
{
// If one of the pattern tests is failed, then the point cannot be
// inside, therefore, if this is a testForInside = true call, the
// result is false. If this is a testForInside = false call, then
// the result is true.
if (surfaceVolumeTests[s][i] != referenceVolumeTypes_[s])
{
insidePoints[i] = false;
inOutSidePoint[i] = !testForInside;
break;
}
}
}
return insidePoints;
return inOutSidePoint;
}
Foam::Field<bool> Foam::conformationSurfaces::wellInside
(
const pointField& samplePts,
const scalarField& testDistSqr
) const
{
return wellInOutSide(samplePts, testDistSqr, true);
}
......@@ -420,9 +437,7 @@ Foam::Field<bool> Foam::conformationSurfaces::wellOutside
const scalarField& testDistSqr
) const
{
notImplemented("Field<bool> Foam::conformationSurfaces::wellOutside");
return Field<bool>(samplePts.size(), true);
return wellInOutSide(samplePts, testDistSqr, false);
}
......
......@@ -176,6 +176,16 @@ public:
//- Check if point is outside surfaces to conform to
bool outside(const point& samplePt) const;
//- Check if point is closer to the surfaces to conform to than
// testDistSqr, in which case return false, otherwise assess in or
// outside and erturn a result depending on the testForInside flag
Field<bool> wellInOutSide
(
const pointField& samplePts,
const scalarField& testDistSqr,
bool testForInside
) const;
//- Check if point is inside surfaces to conform to by at least
// testDistSqr
Field<bool> wellInside
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -2209,38 +2209,45 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
{
vector sampleNearestVec = (sample - nearestPoint);
scalar magSampleNearestVec = mag(sampleNearestVec);
// Nearest to face interior. Use faceNormal to determine side
scalar c = sampleNearestVec & surf.faceNormals()[nearestFaceI];
// If the sample is essentially on the face, do not check for
// it being perpendicular.
if (magSampleNearestVec > 10*SMALL)
{
c /= magSampleNearestVec*mag(surf.faceNormals()[nearestFaceI]);
if (mag(c) < 0.99)
{
FatalErrorIn("triSurfaceTools::surfaceSide")
<< "nearestPoint identified as being on triangle face "
<< "but vector from nearestPoint to sample is not "
<< "perpendicular to the normal." << nl
<< "sample: " << sample << nl
<< "nearestPoint: " << nearestPoint << nl
<< "sample - nearestPoint: " << sample - nearestPoint << nl
<< "normal: " << surf.faceNormals()[nea