diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C index 4eb19041aa..d964ac12cc 100644 --- a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C +++ b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C @@ -262,6 +262,74 @@ bool Foam::checkWedges } +namespace Foam +{ + //- Default transformation behaviour for position + class transformPositionList + { + public: + + //- Transform patch-based field + void operator() + ( + const coupledPolyPatch& cpp, + List& pts + ) const + { + // Each element of pts is all the points in the face. Convert into + // lists of size cpp to transform. + + List newPts(pts.size()); + forAll(pts, faceI) + { + newPts[faceI].setSize(pts[faceI].size()); + } + + label index = 0; + while (true) + { + label n = 0; + + // Extract for every face the i'th position + pointField ptsAtIndex(pts.size(), vector::zero); + forAll(cpp, faceI) + { + const pointField& facePts = pts[faceI]; + if (facePts.size() > index) + { + ptsAtIndex[faceI] = facePts[index]; + n++; + } + } + + if (n == 0) + { + break; + } + + // Now ptsAtIndex will have for every face either zero or + // the position of the i'th vertex. Transform. + cpp.transformPosition(ptsAtIndex); + + // Extract back from ptsAtIndex into newPts + forAll(cpp, faceI) + { + pointField& facePts = newPts[faceI]; + if (facePts.size() > index) + { + facePts[index] = ptsAtIndex[faceI]; + } + } + + index++; + } + + pts.transfer(newPts); + } + }; +} + + bool Foam::checkCoupledPoints ( const polyMesh& mesh, @@ -273,164 +341,106 @@ bool Foam::checkCoupledPoints const faceList& fcs = mesh.faces(); const polyBoundaryMesh& patches = mesh.boundaryMesh(); - // Check size of faces - label maxSize = 0; - { - labelList nbrSize(fcs.size()-mesh.nInternalFaces(), 0); + // Zero'th point on coupled faces + //pointField nbrZeroPoint(fcs.size()-mesh.nInternalFaces(), vector::max); + List nbrPoints(fcs.size() - mesh.nInternalFaces()); - // Exchange size - forAll(patches, patchI) + // Exchange zero point + forAll(patches, patchI) + { + if (patches[patchI].coupled()) { - if (patches[patchI].coupled()) + const coupledPolyPatch& cpp = refCast + ( + patches[patchI] + ); + + forAll(cpp, i) { - const coupledPolyPatch& cpp = refCast + label bFaceI = cpp.start() + i - mesh.nInternalFaces(); + const face& f = cpp[i]; + nbrPoints[bFaceI].setSize(f.size()); + forAll(f, fp) + { + const point& p0 = p[f[fp]]; + nbrPoints[bFaceI][fp] = p0; + } + } + } + } + syncTools::syncBoundaryFaceList + ( + mesh, + nbrPoints, + eqOp(), + transformPositionList() + ); + + // Compare to local ones. Use same tolerance as for matching + label nErrorFaces = 0; + scalar avgMismatch = 0; + label nCoupledPoints = 0; + + forAll(patches, patchI) + { + if (patches[patchI].coupled()) + { + const coupledPolyPatch& cpp = refCast + ( + patches[patchI] + ); + + if (cpp.owner()) + { + scalarField smallDist ( - patches[patchI] + cpp.calcFaceTol + ( + //cpp.matchTolerance(), + cpp, + cpp.points(), + cpp.faceCentres() + ) ); forAll(cpp, i) { - label bFaceI = cpp.start()+i-mesh.nInternalFaces(); - nbrSize[bFaceI] = cpp[i].size(); - maxSize = max(maxSize, cpp[i].size()); - } - } - } - syncTools::swapBoundaryFaceList(mesh, nbrSize); + label bFaceI = cpp.start() + i - mesh.nInternalFaces(); + const face& f = cpp[i]; - - // Check on owner - label nErrorFaces = 0; - forAll(patches, patchI) - { - if (patches[patchI].coupled()) - { - const coupledPolyPatch& cpp = refCast - ( - patches[patchI] - ); - - if (cpp.owner()) - { - forAll(cpp, i) + if (f.size() != nbrPoints[bFaceI].size()) { - label bFaceI = cpp.start()+i-mesh.nInternalFaces(); + FatalErrorIn + ( + "Foam::checkCoupledPoints\n" + "(\n" + " const polyMesh&, const bool, labelHashSet*\n" + ")\n" + ) << "Local face size : " << f.size() + << " does not equal neighbour face size : " + << nbrPoints[bFaceI].size() + << abort(FatalError); + } - if (cpp[i].size() != nbrSize[bFaceI]) + label fp = 0; + forAll(f, j) + { + const point& p0 = p[f[fp]]; + scalar d = mag(p0 - nbrPoints[bFaceI][j]); + + if (d > smallDist[i]) { if (setPtr) { setPtr->insert(cpp.start()+i); } nErrorFaces++; + break; } - } - } - } - } + avgMismatch += d; + nCoupledPoints++; - reduce(nErrorFaces, sumOp