diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C index 960e35ddd4..8dd2eaf423 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C +++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C @@ -1172,6 +1172,8 @@ Foam::tmp Foam::polyMesh::movePoints ) ).movePoints(points_); } + meshSearchMeshObject::Delete(*this); + meshSearchFACECENTRETETSMeshObject::Delete(*this); return sweptVols; } diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C b/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C index bba7367e8b..4dc64d17d8 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C @@ -95,6 +95,8 @@ void Foam::polyMesh::updateMesh(const mapPolyMesh& mpm) ) ).updateMesh(mpm); } + meshSearchMeshObject::Delete(*this); + meshSearchFACECENTRETETSMeshObject::Delete(*this); } diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files index d8687b042b..2be634334c 100644 --- a/src/meshTools/Make/files +++ b/src/meshTools/Make/files @@ -27,6 +27,8 @@ polyMeshZipUpCells/polyMeshZipUpCells.C primitiveMeshGeometry/primitiveMeshGeometry.C meshSearch/meshSearch.C +meshSearch/meshSearchFACECENTRETETSMeshObject.C +meshSearch/meshSearchMeshObject.C meshTools/meshTools.C diff --git a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C index 2b62b29e34..9950b9b061 100644 --- a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C +++ b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C @@ -26,7 +26,7 @@ License #include "mappedPatchBase.H" #include "addToRunTimeSelectionTable.H" #include "ListListOps.H" -#include "meshSearch.H" +#include "meshSearchMeshObject.H" #include "meshTools.H" #include "OFstream.H" #include "Random.H" @@ -37,6 +37,7 @@ License #include "Time.H" #include "mapDistribute.H" #include "SubField.H" +#include "triPointRef.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -80,64 +81,144 @@ const Foam::NamedEnum // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // +Foam::tmp Foam::mappedPatchBase::facePoints +( + const polyPatch& pp +) const +{ + const polyMesh& mesh = pp.boundaryMesh().mesh(); + const labelList& tetBasePts = mesh.tetBasePtIs(); + const pointField& p = pp.points(); + const pointField& cellCentres = mesh.cellCentres(); + const labelList& faceCells = pp.faceCells(); + + // Initialise to face-centre + tmp tfacePoints(new pointField(patch_.faceCentres())); + pointField& facePoints = tfacePoints(); + + forAll(pp, faceI) + { + const face& f = pp[faceI]; + + if (f.size() <= 3) + { + // Point already on tets of cell. + } + else + { + // Find intersection of (face-centre-decomposition) centre to + // cell-centre with face-diagonal-decomposition triangles. + const point& cc = cellCentres[faceCells[faceI]]; + vector d = facePoints[faceI]-cc; + + const label fp0 = tetBasePts[faceI]; + const point& basePoint = p[f[fp0]]; + + //bool foundPoint = false; + + label fp = f.fcIndex(fp0); + for (label i = 2; i < f.size(); i++) + { + const point& thisPoint = p[f[fp]]; + label nextFp = f.fcIndex(fp); + const point& nextPoint = p[f[nextFp]]; + + const triPointRef tri(basePoint, thisPoint, nextPoint); + pointHit hitInfo = tri.intersection + ( + cc, + d, + intersection::HALF_RAY + ); + + if (hitInfo.hit() && hitInfo.distance() > 0) + { + facePoints[faceI] = hitInfo.hitPoint(); + //foundPoint = true; + break; + } + + fp = nextFp; + } + + //if (!foundPoint) + //{ + // Pout<< "cell:" << faceCells[faceI] << " at:" << cc + // << " face-centre:" << patch_.faceCentres()[faceI] + // << " did not find any intersection." << endl; + //} + } + } + + return tfacePoints; +} + + void Foam::mappedPatchBase::collectSamples ( + const pointField& facePoints, pointField& samples, labelList& patchFaceProcs, labelList& patchFaces, pointField& patchFc ) const { - // Collect all sample points and the faces they come from. - List globalFc(Pstream::nProcs()); - List globalSamples(Pstream::nProcs()); - labelListList globalFaces(Pstream::nProcs()); + { + List globalFc(Pstream::nProcs()); + globalFc[Pstream::myProcNo()] = facePoints; + Pstream::gatherList(globalFc); + Pstream::scatterList(globalFc); + // Rework into straight list + patchFc = ListListOps::combine + ( + globalFc, + accessOp() + ); + } - globalFc[Pstream::myProcNo()] = patch_.faceCentres(); - globalSamples[Pstream::myProcNo()] = samplePoints(); - globalFaces[Pstream::myProcNo()] = identity(patch_.size()); + { + List globalSamples(Pstream::nProcs()); + globalSamples[Pstream::myProcNo()] = samplePoints(facePoints); + Pstream::gatherList(globalSamples); + Pstream::scatterList(globalSamples); + // Rework into straight list + samples = ListListOps::combine + ( + globalSamples, + accessOp() + ); + } - // Distribute to all processors - Pstream::gatherList(globalSamples); - Pstream::scatterList(globalSamples); - Pstream::gatherList(globalFaces); - Pstream::scatterList(globalFaces); - Pstream::gatherList(globalFc); - Pstream::scatterList(globalFc); + { + labelListList globalFaces(Pstream::nProcs()); + globalFaces[Pstream::myProcNo()] = identity(patch_.size()); + // Distribute to all processors + Pstream::gatherList(globalFaces); + Pstream::scatterList(globalFaces); - // Rework into straight list - samples = ListListOps::combine - ( - globalSamples, - accessOp() - ); - patchFaces = ListListOps::combine - ( - globalFaces, - accessOp() - ); - patchFc = ListListOps::combine - ( - globalFc, - accessOp() - ); - - patchFaceProcs.setSize(patchFaces.size()); - labelList nPerProc - ( - ListListOps::subSizes + patchFaces = ListListOps::combine ( globalFaces, accessOp() - ) - ); - label sampleI = 0; - forAll(nPerProc, procI) + ); + } + { - for (label i = 0; i < nPerProc[procI]; i++) + labelList nPerProc(Pstream::nProcs()); + nPerProc[Pstream::myProcNo()] = patch_.size(); + Pstream::gatherList(nPerProc); + Pstream::scatterList(nPerProc); + + patchFaceProcs.setSize(patchFaces.size()); + + label sampleI = 0; + forAll(nPerProc, procI) { - patchFaceProcs[sampleI++] = procI; + for (label i = 0; i < nPerProc[procI]; i++) + { + patchFaceProcs[sampleI++] = procI; + } } } } @@ -173,8 +254,9 @@ void Foam::mappedPatchBase::findSamples << sampleModeNames_[mode_] << " mode." << exit(FatalError); } - // Octree based search engine - meshSearch meshSearchEngine(mesh); + //- Note: face-diagonal decomposition + const meshSearchMeshObject& meshSearchEngine = + meshSearchMeshObject::New(mesh); forAll(samples, sampleI) { @@ -290,8 +372,9 @@ void Foam::mappedPatchBase::findSamples << sampleModeNames_[mode_] << " mode." << exit(FatalError); } - // Octree based search engine - meshSearch meshSearchEngine(mesh); + //- Note: face-diagonal decomposition + const meshSearchMeshObject& meshSearchEngine = + meshSearchMeshObject::New(mesh); forAll(samples, sampleI) { @@ -355,23 +438,6 @@ void Foam::mappedPatchBase::findSamples } } - // Check for samples not being found - forAll(nearest, sampleI) - { - if (!nearest[sampleI].first().hit()) - { - FatalErrorIn - ( - "mappedPatchBase::findSamples" - "(const pointField&, labelList&" - ", labelList&, pointField&)" - ) << "Did not find sample " << samples[sampleI] - << " on any processor of region " << sampleRegion_ - << exit(FatalError); - } - } - - // Convert back into proc+local index sampleProcs.setSize(samples.size()); sampleIndices.setSize(samples.size()); @@ -379,21 +445,40 @@ void Foam::mappedPatchBase::findSamples forAll(nearest, sampleI) { - sampleProcs[sampleI] = nearest[sampleI].second().second(); - sampleIndices[sampleI] = nearest[sampleI].first().index(); - sampleLocations[sampleI] = nearest[sampleI].first().hitPoint(); + if (!nearest[sampleI].first().hit()) + { + sampleProcs[sampleI] = -1; + sampleIndices[sampleI] = -1; + sampleLocations[sampleI] = vector::max; + } + else + { + sampleProcs[sampleI] = nearest[sampleI].second().second(); + sampleIndices[sampleI] = nearest[sampleI].first().index(); + sampleLocations[sampleI] = nearest[sampleI].first().hitPoint(); + } } } void Foam::mappedPatchBase::calcMapping() const { + static bool hasWarned = false; + if (mapPtr_.valid()) { FatalErrorIn("mappedPatchBase::calcMapping() const") << "Mapping already calculated" << exit(FatalError); } + // Get points on face (since cannot use face-centres - might be off + // face-diagonal decomposed tets. + tmp patchPoints(facePoints(patch_)); + + // Get offsetted points + const pointField offsettedPoints = samplePoints(patchPoints()); + + // Do a sanity check // Am I sampling my own patch? This only makes sense for a non-zero // offset. @@ -405,8 +490,10 @@ void Foam::mappedPatchBase::calcMapping() const ); // Check offset - vectorField d(samplePoints()-patch_.faceCentres()); - if (sampleMyself && gAverage(mag(d)) <= ROOTVSMALL) + vectorField d(offsettedPoints-patchPoints()); + bool coincident = (gAverage(mag(d)) <= ROOTVSMALL); + + if (sampleMyself && coincident) { WarningIn ( @@ -438,7 +525,15 @@ void Foam::mappedPatchBase::calcMapping() const labelList patchFaceProcs; labelList patchFaces; pointField patchFc; - collectSamples(samples, patchFaceProcs, patchFaces, patchFc); + collectSamples + ( + patchPoints, + samples, + patchFaceProcs, + patchFaces, + patchFc + ); + // Find processor and cell/face samples are in and actual location. labelList sampleProcs; @@ -446,6 +541,75 @@ void Foam::mappedPatchBase::calcMapping() const pointField sampleLocations; findSamples(samples, sampleProcs, sampleIndices, sampleLocations); + // Check for samples that were not found. This will only happen for + // NEARESTCELL since finds cell containing a location + if (mode_ == NEARESTCELL) + { + label nNotFound = 0; + forAll(sampleProcs, sampleI) + { + if (sampleProcs[sampleI] == -1) + { + nNotFound++; + } + } + reduce(nNotFound, sumOp