Commit 4cdb4a2c authored by laurence's avatar laurence
Browse files

ENH: Ensure checkCoupledPoints works in parallel

parent f941955f
......@@ -262,175 +262,185 @@ bool Foam::checkWedges
}
bool Foam::checkCoupledPoints
(
const polyMesh& mesh,
const bool report,
labelHashSet* setPtr
)
namespace Foam
{
const pointField& p = mesh.points();
const faceList& fcs = mesh.faces();
const polyBoundaryMesh& patches = mesh.boundaryMesh();
// Check size of faces
label maxSize = 0;
//- Default transformation behaviour for position
class transformPositionList
{
labelList nbrSize(fcs.size()-mesh.nInternalFaces(), 0);
public:
// Exchange size
forAll(patches, patchI)
//- Transform patch-based field
void operator()
(
const coupledPolyPatch& cpp,
List<pointField>& pts
) const
{
if (patches[patchI].coupled())
// Each element of pts is all the points in the face. Convert into
// lists of size cpp to transform.
List<pointField> newPts(pts.size());
forAll(pts, faceI)
{
const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
(
patches[patchI]
);
newPts[faceI].setSize(pts[faceI].size());
}
forAll(cpp, i)
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)
{
label bFaceI = cpp.start()+i-mesh.nInternalFaces();
nbrSize[bFaceI] = cpp[i].size();
maxSize = max(maxSize, cpp[i].size());
const pointField& facePts = pts[faceI];
if (facePts.size() > index)
{
ptsAtIndex[faceI] = facePts[index];
n++;
}
}
}
}
syncTools::swapBoundaryFaceList(mesh, nbrSize);
if (n == 0)
{
break;
}
// Check on owner
label nErrorFaces = 0;
forAll(patches, patchI)
{
if (patches[patchI].coupled())
{
const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
(
patches[patchI]
);
// Now ptsAtIndex will have for every face either zero or
// the position of the i'th vertex. Transform.
cpp.transformPosition(ptsAtIndex);
if (cpp.owner())
// Extract back from ptsAtIndex into newPts
forAll(cpp, faceI)
{
forAll(cpp, i)
pointField& facePts = newPts[faceI];
if (facePts.size() > index)
{
label bFaceI = cpp.start()+i-mesh.nInternalFaces();
if (cpp[i].size() != nbrSize[bFaceI])
{
if (setPtr)
{
setPtr->insert(cpp.start()+i);
}
nErrorFaces++;
}
facePts[index] = ptsAtIndex[faceI];
}
}
index++;
}
pts.transfer(newPts);
}
};
}
reduce(nErrorFaces, sumOp<label>());
if (nErrorFaces > 0)
bool Foam::checkCoupledPoints
(
const polyMesh& mesh,
const bool report,
labelHashSet* setPtr
)
{
const pointField& p = mesh.points();
const faceList& fcs = mesh.faces();
const polyBoundaryMesh& patches = mesh.boundaryMesh();
// Zero'th point on coupled faces
//pointField nbrZeroPoint(fcs.size()-mesh.nInternalFaces(), vector::max);
List<pointField> nbrPoints(fcs.size() - mesh.nInternalFaces());
// Exchange zero point
forAll(patches, patchI)
{
if (patches[patchI].coupled())
{
if (report)
const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
(
patches[patchI]
);
forAll(cpp, i)
{
Info<< " **Error in coupled faces: "
<< nErrorFaces
<< " faces have different size "
<< " compared to their coupled equivalent." << endl;
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;
}
}
return true;
}
reduce(maxSize, maxOp<label>());
}
syncTools::syncBoundaryFaceList
(
mesh,
nbrPoints,
eqOp<pointField>(),
transformPositionList()
);
// Compare to local ones. Use same tolerance as for matching
label nErrorFaces = 0;
scalar avgMismatch = 0;
label nCoupledPoints = 0;
for (label index = 0; index < maxSize; index++)
forAll(patches, patchI)
{
// point at index on coupled faces
pointField nbrPoint(fcs.size()-mesh.nInternalFaces(), vector::max);
// Exchange point
forAll(patches, patchI)
if (patches[patchI].coupled())
{
if (patches[patchI].coupled())
const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
(
patches[patchI]
);
if (cpp.owner())
{
const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
scalarField smallDist
(
patches[patchI]
cpp.calcFaceTol
(
//cpp.matchTolerance(),
cpp,
cpp.points(),
cpp.faceCentres()
)
);
forAll(cpp, i)
{
label bFaceI = cpp.start() + i - mesh.nInternalFaces();
const face& f = cpp[i];
if (f.size() > index)
if (f.size() != nbrPoints[bFaceI].size())
{
label bFaceI = cpp.start()+i-mesh.nInternalFaces();
nbrPoint[bFaceI] = p[f[index]];
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);
}
}
}
}
syncTools::swapBoundaryFacePositions(mesh, nbrPoint);
// Compare to local ones. Use same tolerance as for matching
forAll(patches, patchI)
{
if (patches[patchI].coupled())
{
const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
(
patches[patchI]
);
if (cpp.owner())
{
scalarField smallDist
(
cpp.calcFaceTol
(
//cpp.matchTolerance(),
cpp,
cpp.points(),
cpp.faceCentres()
)
);
forAll(cpp, i)
label fp = 0;
forAll(f, j)
{
const face& f = cpp[i];
if (f.size() > index)
{
label bFaceI = cpp.start()+i-mesh.nInternalFaces();
label reverseIndex = (f.size()-index)%f.size();
scalar d = mag(p[f[reverseIndex]]-nbrPoint[bFaceI]);
const point& p0 = p[f[fp]];
scalar d = mag(p0 - nbrPoints[bFaceI][j]);
if (d > smallDist[i])
if (d > smallDist[i])
{
if (setPtr)
{
if (setPtr)
{
// Avoid duplicate counting of faces
if (setPtr->insert(cpp.start()+i))
{
nErrorFaces++;
}
}
else
{
// No checking on duplicates
nErrorFaces++;
}
setPtr->insert(cpp.start()+i);
}
avgMismatch += d;
nCoupledPoints++;
nErrorFaces++;
break;
}
avgMismatch += d;
nCoupledPoints++;
fp = f.rcIndex(fp);
}
}
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment