diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolume.C b/src/meshTools/tetOverlapVolume/tetOverlapVolume.C index c430a4803178313881b06c82845abf9e89cbd0d1..c0f192f7de4a6b88a3e25708d2c5155452b6e2a3 100644 --- a/src/meshTools/tetOverlapVolume/tetOverlapVolume.C +++ b/src/meshTools/tetOverlapVolume/tetOverlapVolume.C @@ -130,6 +130,135 @@ Foam::treeBoundBox Foam::tetOverlapVolume::pyrBb // * * * * * * * * * * * Public Member Functions * * * * * * * * * * * * * // +bool Foam::tetOverlapVolume::cellCellOverlapMinDecomp +( + const primitiveMesh& meshA, + const label cellAI, + const primitiveMesh& meshB, + const label cellBI, + const treeBoundBox& cellBbB, + const scalar threshold +) const +{ + const cell& cFacesA = meshA.cells()[cellAI]; + const point& ccA = meshA.cellCentres()[cellAI]; + + const cell& cFacesB = meshB.cells()[cellBI]; + const point& ccB = meshB.cellCentres()[cellBI]; + + scalar vol = 0.0; + + forAll(cFacesA, cFA) + { + label faceAI = cFacesA[cFA]; + + const face& fA = meshA.faces()[faceAI]; + const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA); + if (!pyrA.overlaps(cellBbB)) + { + continue; + } + + bool ownA = (meshA.faceOwner()[faceAI] == cellAI); + + label tetBasePtAI = 0; + + const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]]; + + for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++) + { + label facePtAI = (tetPtI + tetBasePtAI) % fA.size(); + label otherFacePtAI = fA.fcIndex(facePtAI); + + label pt0I = -1; + label pt1I = -1; + + if (ownA) + { + pt0I = fA[facePtAI]; + pt1I = fA[otherFacePtAI]; + } + else + { + pt0I = fA[otherFacePtAI]; + pt1I = fA[facePtAI]; + } + + const tetPoints tetA + ( + ccA, + tetBasePtA, + meshA.points()[pt0I], + meshA.points()[pt1I] + ); + const treeBoundBox tetABb(tetA.bounds()); + + + // Loop over tets of cellB + forAll(cFacesB, cFB) + { + label faceBI = cFacesB[cFB]; + + const face& fB = meshB.faces()[faceBI]; + const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB); + if (!pyrB.overlaps(pyrA)) + { + continue; + } + + bool ownB = (meshB.faceOwner()[faceBI] == cellBI); + + label tetBasePtBI = 0; + + const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]]; + + for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++) + { + label facePtBI = (tetPtI + tetBasePtBI) % fB.size(); + label otherFacePtBI = fB.fcIndex(facePtBI); + + label pt0I = -1; + label pt1I = -1; + + if (ownB) + { + pt0I = fB[facePtBI]; + pt1I = fB[otherFacePtBI]; + } + else + { + pt0I = fB[otherFacePtBI]; + pt1I = fB[facePtBI]; + } + + const tetPoints tetB + ( + ccB, + tetBasePtB, + meshB.points()[pt0I], + meshB.points()[pt1I] + ); + + if (!tetB.bounds().overlaps(tetABb)) + { + continue; + } + + vol += tetTetOverlapVol(tetA, tetB); + + if (vol > threshold) + { + return true; + } + } + } + } + } + + return false; +} + + Foam::scalar Foam::tetOverlapVolume::cellCellOverlapVolumeMinDecomp ( const primitiveMesh& meshA, diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolume.H b/src/meshTools/tetOverlapVolume/tetOverlapVolume.H index 631ed9a50d93ec139611eb0e7d02a51002ec84cf..58646f2e5a7ce3693391bbdb6d01acfb6f2284c3 100644 --- a/src/meshTools/tetOverlapVolume/tetOverlapVolume.H +++ b/src/meshTools/tetOverlapVolume/tetOverlapVolume.H @@ -48,7 +48,7 @@ class polyMesh; class tetPoints; /*---------------------------------------------------------------------------*\ - Class tetOverlapVolume Declaration + Class tetOverlapVolume Declaration \*---------------------------------------------------------------------------*/ class tetOverlapVolume @@ -94,6 +94,16 @@ public: const label cellBI ) const; + //- Return true if olverlap volume is greater than threshold + bool cellCellOverlapMinDecomp + ( + const primitiveMesh& meshA, + const label cellAI, + const primitiveMesh& meshB, + const label cellBI, + const treeBoundBox& cellBbB, + const scalar threshold = 0.0 + ) const; //- Calculates the overlap volume scalar cellCellOverlapVolumeMinDecomp