diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C index 8bec6c4ae8fd57276ec5513d1e0d5b4c5040c4cb..9f466dd13bee979628b097588ff71512eddea9dd 100644 --- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C +++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C @@ -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." diff --git a/src/mesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.C b/src/mesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.C index 19169d04f572233cc8127292536b1f8588a81c7f..bda7d08e70f3b7e197f3b63179be5070cbe39ed0 100644 --- a/src/mesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.C +++ b/src/mesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.C @@ -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 diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H index c8647b4bd316a6966c5a3c177f0ee3fac646e7ae..1104440dddefddb32a115dfc494ef104d4adc789 100644 --- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H +++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H @@ -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; }; diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C index 2068ba3e10dc05a3202dc39ee518bff4d150a7e5..0371b6adca257c7df59a44ce21f1ec8d075b592f 100644 --- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C +++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C @@ -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); diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C index bd2d99932eb9c6146708d2737a5cd81629d63f67..fefc5546d6b8974cc5c05f41f602ce01ada001b1 100644 --- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C +++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C @@ -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(); + } +} + + // ************************************************************************* // diff --git a/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C b/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C index fe71457794e236ee7c509b436caeaae7432327b5..527a15023ca582e2e5df246c9116ebcf631d765f 100644 --- a/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C +++ b/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C @@ -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); } diff --git a/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H b/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H index a7a952eeddd90aeb05499f27abd432eadaafca45..ad472812b3a9191fddaffd5f80bdb8d51446f6ad 100644 --- a/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H +++ b/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H @@ -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 diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C index cf5d76cb95f1c351deef876d467512a0f17f9ced..d5819eb2f7e6f8c6f6618d5716f00b28232fe843 100644 --- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C +++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C @@ -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()[nearestFaceI] << nl - << "mag(sample - nearestPoint): " - << mag(sample - nearestPoint) << nl - << "normalised dot product: " << c << nl - << "triangle vertices: " << nl - << " " << points[f[0]] << nl - << " " << points[f[1]] << nl - << " " << points[f[2]] << nl - << abort(FatalError); - } - } + // if (debug) + // { + // scalar magSampleNearestVec = mag(sampleNearestVec); + + // if (magSampleNearestVec > SMALL) + // { + // c /= + // magSampleNearestVec + // *mag(surf.faceNormals()[nearestFaceI]); + + // if (mag(c) < 0.99) + // { + // WarningIn("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()[nearestFaceI] << nl + // << "mag(sample - nearestPoint): " + // << mag(sample - nearestPoint) << nl + // << "normalised dot product: " << c << nl + // << "triangle vertices: " << nl + // << " " << points[f[0]] << nl + // << " " << points[f[1]] << nl + // << " " << points[f[2]] << nl + // << endl;; + // } + // } + // } if (c > 0) {