diff --git a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C index ef0ee2d85b490af5e5a6d6d42f905f6d6e8091c1..b2337c67bf8abf4c77ac97a5e4cbb57364fe3f34 100644 --- a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C +++ b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C @@ -259,9 +259,9 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict) ); wallDistCalc.iterate(nMedialAxisIter); - label nUnvisit = returnReduce + const label nUnvisit = returnReduce ( - wallDistCalc.getUnsetPoints(), + wallDistCalc.nUnvisitedPoints(), sumOp<label>() ); diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C index dc9d2ba85ad3e2612e46353ff4fa25a08fb90d26..28271824c7cd44783de7a1374d07f5902884231c 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C +++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C @@ -108,6 +108,16 @@ Foam::meshRefinement::writeType Foam::meshRefinement::writeLevel_; Foam::meshRefinement::outputType Foam::meshRefinement::outputLevel_; +// Inside/outside test for polyMesh:.findCell() +// 2.4.x : default = polyMesh::FACE_DIAG_TRIS +// 1712 : default = polyMesh::CELL_TETS +// +// - CELL_TETS is better with concave cells, but much slower. +// - use faster method (FACE_DIAG_TRIS) here + +static const Foam::polyMesh::cellDecomposition + findCellMode(Foam::polyMesh::FACE_DIAG_TRIS); + // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // @@ -2198,18 +2208,7 @@ Foam::label Foam::meshRefinement::findRegion // Force calculation of base points (needs to be synchronised) (void)mesh.tetBasePtIs(); - label celli = mesh.findCell(p); - //if (celli != -1) - //{ - // Pout<< "findRegion : Found point:" << p << " in cell " << celli - // << " at:" << mesh.cellCentres()[celli] << endl; - //} - //else - //{ - // Pout<< "findRegion : Found point:" << p << " in cell " << celli - // << endl; - //} - + label celli = mesh.findCell(p, findCellMode); if (celli != -1) { regioni = cellToRegion[celli]; @@ -2219,7 +2218,7 @@ Foam::label Foam::meshRefinement::findRegion if (regioni == -1) { // See if we can perturb a bit - celli = mesh.findCell(p+perturbVec); + celli = mesh.findCell(p+perturbVec, findCellMode); if (celli != -1) { regioni = cellToRegion[celli]; @@ -2228,61 +2227,12 @@ Foam::label Foam::meshRefinement::findRegion } return regioni; } -//XXXXXXXX -//Foam::labelList Foam::meshRefinement::findRegion -//( -// const polyMesh& mesh, -// const labelList& cellToRegion, -// const vector& perturbVec, -// const pointField& pts -//) -//{ -// labelList regions(pts.size(), -1); -// -// forAll(pts, i) -// { -// label celli = mesh.findCell(pts[i]); -// if (celli != -1) -// { -// regions[i] = cellToRegion[celli]; -// } -// reduce(regions[i], maxOp<label>()); -// -// if (regions[i] == -1) -// { -// // See if we can perturb a bit -// celli = mesh.findCell(pts[i]+perturbVec); -// if (celli != -1) -// { -// regions[i] = cellToRegion[celli]; -// } -// reduce(regions[i], maxOp<label>()); -// } -// } -// -// forAll(regions, i) -// { -// if (regions[i] == -1) -// { -// FatalErrorInFunction -// << "Point " << pts[i] -// << " is not inside the mesh." << nl -// << "Bounding box of the mesh:" << mesh.bounds() -// //<< "All points " << pts -// //<< " with corresponding regions " << regions -// << exit(FatalError); -// } -// } -// -// return regions; -//} -//XXXXXXXX // Modify cellRegion to be consistent with locationsInMesh. // - all regions not in locationsInMesh are set to -1 // - check that all regions inside locationsOutsideMesh are set to -1 -void Foam::meshRefinement::findRegions +Foam::label Foam::meshRefinement::findRegions ( const polyMesh& mesh, const vector& perturbVec, @@ -2352,14 +2302,23 @@ void Foam::meshRefinement::findRegions } + label nRemove = 0; + // Now update cellRegion to -1 for unreachable cells forAll(insideCell, celli) { - if (!insideCell[celli]) + if (!insideCell.test(celli)) { cellRegion[celli] = -1; + ++nRemove; + } + else if (cellRegion[celli] == -1) + { + ++nRemove; } } + + return nRemove; } @@ -2382,7 +2341,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMeshRegions regionSplit cellRegion(mesh_, blockedFace); - findRegions + label nRemove = findRegions ( mesh_, mergeDistance_*vector(1,1,1), // perturbVec @@ -2396,7 +2355,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMeshRegions // ~~~~~~ // Get cells to remove - DynamicList<label> cellsToRemove(mesh_.nCells()); + DynamicList<label> cellsToRemove(nRemove); forAll(cellRegion, celli) { if (cellRegion[celli] == -1) diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H index 3c5c03f74cfa60212f83fe4878f4f72de0deef1f..fd36951239a6e39718203944af3400b0e00b5f91 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H +++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H @@ -1226,7 +1226,9 @@ public: const point& p ); - static void findRegions + //- Find regions points are in. + // \return number of cells to be removed + static label findRegions ( const polyMesh&, const vector& perturbVec, diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C index 50142bbbc92c13653f1e0d219df1e139e0fe9f7c..cda200b52b362bff990212ab54c4f6df00503fa4 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C +++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C @@ -1570,7 +1570,7 @@ Foam::label Foam::meshRefinement::markSmallFeatureRefinement ); //- Option 1: use octree nearest searching inside polyMesh - //label cellI = mesh_.findCell(pt); + //label cellI = mesh_.findCell(pt, polyMesh::CELL_TETS); //- Option 2: use octree 'inside' searching inside polyMesh. Is // much faster. diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.C index c4a7ed6c1328faad6830dfb75f3a0e48282f6d33..836f8ef826511be152b57e772b488a99f00a44ee 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.C +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.C @@ -191,7 +191,7 @@ Foam::labelList Foam::refinementParameters::findCells { const point& location = locations[i]; - label localCellI = mesh.findCell(location); + label localCellI = mesh.findCell(location, polyMesh::FACE_DIAG_TRIS); label globalCellI = -1; diff --git a/src/meshTools/algorithms/MeshWave/FaceCellWave.C b/src/meshTools/algorithms/MeshWave/FaceCellWave.C index 8fcaa8cbe6b876f056616a64cdd16250e8e1360a..67226d5d51f5466df91bb2857c67b237b1a23e49 100644 --- a/src/meshTools/algorithms/MeshWave/FaceCellWave.C +++ b/src/meshTools/algorithms/MeshWave/FaceCellWave.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -122,11 +122,11 @@ bool Foam::FaceCellWave<Type, TrackingData>::updateCell // - changedCell_, changedCells_ // - statistics: nEvals_, nUnvisitedCells_ - nEvals_++; + ++nEvals_; - bool wasValid = cellInfo.valid(td_); + const bool wasValid = cellInfo.valid(td_); - bool propagate = + const bool propagate = cellInfo.updateCell ( mesh_, @@ -139,9 +139,8 @@ bool Foam::FaceCellWave<Type, TrackingData>::updateCell if (propagate) { - if (!changedCell_[celli]) + if (changedCell_.set(celli)) { - changedCell_.set(celli); changedCells_.append(celli); } } @@ -171,11 +170,11 @@ bool Foam::FaceCellWave<Type, TrackingData>::updateFace // - changedFace_, changedFaces_, // - statistics: nEvals_, nUnvisitedFaces_ - nEvals_++; + ++nEvals_; - bool wasValid = faceInfo.valid(td_); + const bool wasValid = faceInfo.valid(td_); - bool propagate = + const bool propagate = faceInfo.updateFace ( mesh_, @@ -218,11 +217,11 @@ bool Foam::FaceCellWave<Type, TrackingData>::updateFace // - changedFace_, changedFaces_, // - statistics: nEvals_, nUnvisitedFaces_ - nEvals_++; + ++nEvals_; - bool wasValid = faceInfo.valid(td_); + const bool wasValid = faceInfo.valid(td_); - bool propagate = + const bool propagate = faceInfo.updateFace ( mesh_, @@ -262,8 +261,8 @@ void Foam::FaceCellWave<Type, TrackingData>::checkCyclic forAll(patch, patchFacei) { - label i1 = patch.start() + patchFacei; - label i2 = nbrPatch.start() + patchFacei; + const label i1 = patch.start() + patchFacei; + const label i2 = nbrPatch.start() + patchFacei; if ( @@ -299,9 +298,9 @@ template<class Type, class TrackingData> template<class PatchType> bool Foam::FaceCellWave<Type, TrackingData>::hasPatch() const { - forAll(mesh_.boundaryMesh(), patchi) + for (const polyPatch& p : mesh_.boundaryMesh()) { - if (isA<PatchType>(mesh_.boundaryMesh()[patchi])) + if (isA<PatchType>(p)) { return true; } @@ -313,15 +312,39 @@ bool Foam::FaceCellWave<Type, TrackingData>::hasPatch() const template<class Type, class TrackingData> void Foam::FaceCellWave<Type, TrackingData>::setFaceInfo ( - const labelList& changedFaces, + const label facei, + const Type& faceInfo +) +{ + const bool wasValid = allFaceInfo_[facei].valid(td_); + + // Copy info for facei + allFaceInfo_[facei] = faceInfo; + + // Maintain count of unset faces + if (!wasValid && allFaceInfo_[facei].valid(td_)) + { + --nUnvisitedFaces_; + } + + // Mark facei as visited and changed (both on list and on face itself) + changedFace_.set(facei); + changedFaces_.append(facei); +} + + +template<class Type, class TrackingData> +void Foam::FaceCellWave<Type, TrackingData>::setFaceInfo +( + const labelUList& changedFaces, const List<Type>& changedFacesInfo ) { forAll(changedFaces, changedFacei) { - label facei = changedFaces[changedFacei]; + const label facei = changedFaces[changedFacei]; - bool wasValid = allFaceInfo_[facei].valid(td_); + const bool wasValid = allFaceInfo_[facei].valid(td_); // Copy info for facei allFaceInfo_[facei] = changedFacesInfo[changedFacei]; @@ -333,7 +356,6 @@ void Foam::FaceCellWave<Type, TrackingData>::setFaceInfo } // Mark facei as changed, both on list and on face itself. - changedFace_.set(facei); changedFaces_.append(facei); } @@ -345,29 +367,29 @@ void Foam::FaceCellWave<Type, TrackingData>::mergeFaceInfo ( const polyPatch& patch, const label nFaces, - const labelList& changedFaces, + const labelUList& changedFaces, const List<Type>& changedFacesInfo ) { // Merge face information into member data - for (label changedFacei = 0; changedFacei < nFaces; changedFacei++) + for (label changedFacei = 0; changedFacei < nFaces; ++changedFacei) { - const Type& neighbourWallInfo = changedFacesInfo[changedFacei]; - label patchFacei = changedFaces[changedFacei]; + const Type& newInfo = changedFacesInfo[changedFacei]; + const label patchFacei = changedFaces[changedFacei]; - label meshFacei = patch.start() + patchFacei; + const label meshFacei = patch.start() + patchFacei; - Type& currentWallInfo = allFaceInfo_[meshFacei]; + Type& currInfo = allFaceInfo_[meshFacei]; - if (!currentWallInfo.equal(neighbourWallInfo, td_)) + if (!currInfo.equal(newInfo, td_)) { updateFace ( meshFacei, - neighbourWallInfo, + newInfo, propagationTol_, - currentWallInfo + currInfo ); } } @@ -387,22 +409,22 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::getChangedPatchFaces // Construct compact patchFace change arrays for a (slice of a) single // patch. changedPatchFaces in local patch numbering. // Return length of arrays. - label nChangedPatchFaces = 0; + label nChanged = 0; - for (label i = 0; i < nFaces; i++) + for (label i = 0; i < nFaces; ++i) { - label patchFacei = i + startFacei; - - label meshFacei = patch.start() + patchFacei; + const label patchFacei = i + startFacei; + const label meshFacei = patch.start() + patchFacei; if (changedFace_.test(meshFacei)) { - changedPatchFaces[nChangedPatchFaces] = patchFacei; - changedPatchFacesInfo[nChangedPatchFaces] = allFaceInfo_[meshFacei]; - nChangedPatchFaces++; + changedPatchFaces[nChanged] = patchFacei; + changedPatchFacesInfo[nChanged] = allFaceInfo_[meshFacei]; + ++nChanged; } } - return nChangedPatchFaces; + + return nChanged; } @@ -411,7 +433,7 @@ void Foam::FaceCellWave<Type, TrackingData>::leaveDomain ( const polyPatch& patch, const label nFaces, - const labelList& faceLabels, + const labelUList& faceLabels, List<Type>& faceInfo ) const { @@ -419,11 +441,11 @@ void Foam::FaceCellWave<Type, TrackingData>::leaveDomain const vectorField& fc = mesh_.faceCentres(); - for (label i = 0; i < nFaces; i++) + for (label i = 0; i < nFaces; ++i) { - label patchFacei = faceLabels[i]; + const label patchFacei = faceLabels[i]; + const label meshFacei = patch.start() + patchFacei; - label meshFacei = patch.start() + patchFacei; faceInfo[i].leaveDomain(mesh_, patch, patchFacei, fc[meshFacei], td_); } } @@ -434,7 +456,7 @@ void Foam::FaceCellWave<Type, TrackingData>::enterDomain ( const polyPatch& patch, const label nFaces, - const labelList& faceLabels, + const labelUList& faceLabels, List<Type>& faceInfo ) const { @@ -442,11 +464,11 @@ void Foam::FaceCellWave<Type, TrackingData>::enterDomain const vectorField& fc = mesh_.faceCentres(); - for (label i = 0; i < nFaces; i++) + for (label i = 0; i < nFaces; ++i) { - label patchFacei = faceLabels[i]; + const label patchFacei = faceLabels[i]; + const label meshFacei = patch.start() + patchFacei; - label meshFacei = patch.start() + patchFacei; faceInfo[i].enterDomain(mesh_, patch, patchFacei, fc[meshFacei], td_); } } @@ -466,14 +488,14 @@ void Foam::FaceCellWave<Type, TrackingData>::transform { const tensor& T = rotTensor[0]; - for (label facei = 0; facei < nFaces; facei++) + for (label facei = 0; facei < nFaces; ++facei) { faceInfo[facei].transform(mesh_, T, td_); } } else { - for (label facei = 0; facei < nFaces; facei++) + for (label facei = 0; facei < nFaces; ++facei) { faceInfo[facei].transform(mesh_, rotTensor[facei], td_); } @@ -490,10 +512,10 @@ void Foam::FaceCellWave<Type, TrackingData>::offset labelList& faces ) { - // Offset mesh face. Used for transferring from one cyclic half to the - // other. + // Offset mesh face. + // Used for transferring from one cyclic half to the other. - for (label facei = 0; facei < nFaces; facei++) + for (label facei = 0; facei < nFaces; ++facei) { faces[facei] += cycOffset; } @@ -514,10 +536,8 @@ void Foam::FaceCellWave<Type, TrackingData>::handleProcPatches() PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking); - forAll(procPatches, i) + for (const label patchi : procPatches) { - label patchi = procPatches[i]; - const processorPolyPatch& procPatch = refCast<const processorPolyPatch>(mesh_.boundaryMesh()[patchi]); @@ -564,10 +584,8 @@ void Foam::FaceCellWave<Type, TrackingData>::handleProcPatches() // Receive all - forAll(procPatches, i) + for (const label patchi : procPatches) { - label patchi = procPatches[i]; - const processorPolyPatch& procPatch = refCast<const processorPolyPatch>(mesh_.boundaryMesh()[patchi]); @@ -625,14 +643,14 @@ void Foam::FaceCellWave<Type, TrackingData>::handleCyclicPatches() { // Transfer information across cyclic halves. - forAll(mesh_.boundaryMesh(), patchi) + for (const polyPatch& patch : mesh_.boundaryMesh()) { - const polyPatch& patch = mesh_.boundaryMesh()[patchi]; - if (isA<cyclicPolyPatch>(patch)) { - const cyclicPolyPatch& nbrPatch = - refCast<const cyclicPolyPatch>(patch).neighbPatch(); + const cyclicPolyPatch& cycPatch = + refCast<const cyclicPolyPatch>(patch); + + const cyclicPolyPatch& nbrPatch = cycPatch.neighbPatch(); // Allocate buffers label nReceiveFaces; @@ -658,12 +676,9 @@ void Foam::FaceCellWave<Type, TrackingData>::handleCyclicPatches() receiveFacesInfo ); - const cyclicPolyPatch& cycPatch = - refCast<const cyclicPolyPatch>(patch); - if (!cycPatch.parallel()) { - // received data from other half + // Received data from other half transform ( cycPatch.forwardT(), @@ -674,7 +689,8 @@ void Foam::FaceCellWave<Type, TrackingData>::handleCyclicPatches() if (debug & 2) { - Pout<< " Cyclic patch " << patchi << ' ' << cycPatch.name() + Pout<< " Cyclic patch " + << cycPatch.index() << ' ' << cycPatch.name() << " Changed : " << nReceiveFaces << endl; } @@ -711,21 +727,18 @@ void Foam::FaceCellWave<Type, TrackingData>::handleAMICyclicPatches() { // Transfer information across cyclicAMI halves. - forAll(mesh_.boundaryMesh(), patchi) + for (const polyPatch& patch : mesh_.boundaryMesh()) { - const polyPatch& patch = mesh_.boundaryMesh()[patchi]; - if (isA<cyclicAMIPolyPatch>(patch)) { const cyclicAMIPolyPatch& cycPatch = refCast<const cyclicAMIPolyPatch>(patch); + const cyclicAMIPolyPatch& nbrPatch = cycPatch.neighbPatch(); + List<Type> receiveInfo; { - const cyclicAMIPolyPatch& nbrPatch = - refCast<const cyclicAMIPolyPatch>(patch).neighbPatch(); - // Get nbrPatch data (so not just changed faces) typename List<Type>::subList sendInfo ( @@ -787,22 +800,20 @@ void Foam::FaceCellWave<Type, TrackingData>::handleAMICyclicPatches() // Merge into global storage forAll(receiveInfo, i) { - label meshFacei = cycPatch.start()+i; + const label meshFacei = cycPatch.start()+i; - Type& currentWallInfo = allFaceInfo_[meshFacei]; + const Type& newInfo = receiveInfo[i]; - if - ( - receiveInfo[i].valid(td_) - && !currentWallInfo.equal(receiveInfo[i], td_) - ) + Type& currInfo = allFaceInfo_[meshFacei]; + + if (newInfo.valid(td_) && !currInfo.equal(newInfo, td_)) { updateFace ( meshFacei, - receiveInfo[i], + newInfo, propagationTol_, - currentWallInfo + currInfo ); } } @@ -814,73 +825,50 @@ void Foam::FaceCellWave<Type, TrackingData>::handleAMICyclicPatches() template<class Type, class TrackingData> void Foam::FaceCellWave<Type, TrackingData>::handleExplicitConnections() { - // Collect changed information - - DynamicList<label> f0Baffle(explicitConnections_.size()); - DynamicList<Type> f0Info(explicitConnections_.size()); - - DynamicList<label> f1Baffle(explicitConnections_.size()); - DynamicList<Type> f1Info(explicitConnections_.size()); + changedBaffles_.clear(); - forAll(explicitConnections_, connI) + // Collect all/any changed information touching a baffle + for (const labelPair& baffle : explicitConnections_) { - const labelPair& baffle = explicitConnections_[connI]; + const label f0 = baffle.first(); + const label f1 = baffle.second(); - label f0 = baffle[0]; if (changedFace_.test(f0)) { - f0Baffle.append(connI); - f0Info.append(allFaceInfo_[f0]); + // f0 changed. Update information on f1. + changedBaffles_.append(taggedInfoType(f1, allFaceInfo_[f0])); } - label f1 = baffle[1]; if (changedFace_.test(f1)) { - f1Baffle.append(connI); - f1Info.append(allFaceInfo_[f1]); + // f1 changed. Update information on f0. + changedBaffles_.append(taggedInfoType(f0, allFaceInfo_[f1])); } } // Update other side with changed information - forAll(f1Info, index) + for (const taggedInfoType& updated : changedBaffles_) { - const labelPair& baffle = explicitConnections_[f1Baffle[index]]; + const label tgtFace = updated.first; + const Type& newInfo = updated.second; - label f0 = baffle[0]; - Type& currentWallInfo = allFaceInfo_[f0]; + Type& currInfo = allFaceInfo_[tgtFace]; - if (!currentWallInfo.equal(f1Info[index], td_)) + if (!currInfo.equal(newInfo, td_)) { updateFace ( - f0, - f1Info[index], + tgtFace, + newInfo, propagationTol_, - currentWallInfo + currInfo ); } } - forAll(f0Info, index) - { - const labelPair& baffle = explicitConnections_[f0Baffle[index]]; - - label f1 = baffle[1]; - Type& currentWallInfo = allFaceInfo_[f1]; - - if (!currentWallInfo.equal(f0Info[index], td_)) - { - updateFace - ( - f1, - f0Info[index], - propagationTol_, - currentWallInfo - ); - } - } + changedBaffles_.clear(); } @@ -896,14 +884,15 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave ) : mesh_(mesh), - explicitConnections_(0), + explicitConnections_(), allFaceInfo_(allFaceInfo), allCellInfo_(allCellInfo), td_(td), changedFace_(mesh_.nFaces(), false), - changedFaces_(mesh_.nFaces()), changedCell_(mesh_.nCells(), false), + changedFaces_(mesh_.nFaces()), changedCells_(mesh_.nCells()), + changedBaffles_(2*explicitConnections_.size()), hasCyclicPatches_(hasPatch<cyclicPolyPatch>()), hasCyclicAMIPatches_ ( @@ -920,12 +909,11 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave ) { FatalErrorInFunction - << "face and cell storage not the size of mesh faces, cells:" - << endl - << " allFaceInfo :" << allFaceInfo.size() << endl - << " mesh_.nFaces():" << mesh_.nFaces() << endl - << " allCellInfo :" << allCellInfo.size() << endl - << " mesh_.nCells():" << mesh_.nCells() + << "face and cell storage not the size of mesh faces, cells:" << nl + << " allFaceInfo :" << allFaceInfo.size() << nl + << " mesh_.nFaces():" << mesh_.nFaces() << nl + << " allCellInfo :" << allCellInfo.size() << nl + << " mesh_.nCells():" << mesh_.nCells() << endl << exit(FatalError); } } @@ -935,7 +923,7 @@ template<class Type, class TrackingData> Foam::FaceCellWave<Type, TrackingData>::FaceCellWave ( const polyMesh& mesh, - const labelList& changedFaces, + const labelUList& changedFaces, const List<Type>& changedFacesInfo, UList<Type>& allFaceInfo, UList<Type>& allCellInfo, @@ -944,14 +932,15 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave ) : mesh_(mesh), - explicitConnections_(0), + explicitConnections_(), allFaceInfo_(allFaceInfo), allCellInfo_(allCellInfo), td_(td), changedFace_(mesh_.nFaces(), false), - changedFaces_(mesh_.nFaces()), changedCell_(mesh_.nCells(), false), + changedFaces_(mesh_.nFaces()), changedCells_(mesh_.nCells()), + changedBaffles_(2*explicitConnections_.size()), hasCyclicPatches_(hasPatch<cyclicPolyPatch>()), hasCyclicAMIPatches_ ( @@ -968,12 +957,11 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave ) { FatalErrorInFunction - << "face and cell storage not the size of mesh faces, cells:" - << endl - << " allFaceInfo :" << allFaceInfo.size() << endl - << " mesh_.nFaces():" << mesh_.nFaces() << endl - << " allCellInfo :" << allCellInfo.size() << endl - << " mesh_.nCells():" << mesh_.nCells() + << "face and cell storage not the size of mesh faces, cells:" << nl + << " allFaceInfo :" << allFaceInfo.size() << nl + << " mesh_.nFaces():" << mesh_.nFaces() << nl + << " allCellInfo :" << allCellInfo.size() << nl + << " mesh_.nCells():" << mesh_.nCells() << endl << exit(FatalError); } @@ -981,14 +969,14 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave setFaceInfo(changedFaces, changedFacesInfo); // Iterate until nothing changes - label iter = iterate(maxIter); + const label iter = iterate(maxIter); if ((maxIter > 0) && (iter >= maxIter)) { FatalErrorInFunction - << "Maximum number of iterations reached. Increase maxIter." << endl - << " maxIter:" << maxIter << endl - << " nChangedCells:" << changedCells_.size() << endl + << "Maximum number of iterations reached. Increase maxIter." << nl + << " maxIter:" << maxIter << nl + << " nChangedCells:" << changedCells_.size() << nl << " nChangedFaces:" << changedFaces_.size() << endl << exit(FatalError); } @@ -1001,7 +989,7 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave const polyMesh& mesh, const labelPairList& explicitConnections, const bool handleCyclicAMI, - const labelList& changedFaces, + const labelUList& changedFaces, const List<Type>& changedFacesInfo, UList<Type>& allFaceInfo, UList<Type>& allCellInfo, @@ -1015,9 +1003,10 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave allCellInfo_(allCellInfo), td_(td), changedFace_(mesh_.nFaces(), false), - changedFaces_(mesh_.nFaces()), changedCell_(mesh_.nCells(), false), + changedFaces_(mesh_.nFaces()), changedCells_(mesh_.nCells()), + changedBaffles_(2*explicitConnections_.size()), hasCyclicPatches_(hasPatch<cyclicPolyPatch>()), hasCyclicAMIPatches_ ( @@ -1035,12 +1024,11 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave ) { FatalErrorInFunction - << "face and cell storage not the size of mesh faces, cells:" - << endl - << " allFaceInfo :" << allFaceInfo.size() << endl - << " mesh_.nFaces():" << mesh_.nFaces() << endl - << " allCellInfo :" << allCellInfo.size() << endl - << " mesh_.nCells():" << mesh_.nCells() + << "face and cell storage not the size of mesh faces, cells:" << nl + << " allFaceInfo :" << allFaceInfo.size() << nl + << " mesh_.nFaces():" << mesh_.nFaces() << nl + << " allCellInfo :" << allCellInfo.size() << nl + << " mesh_.nCells():" << mesh_.nCells() << endl << exit(FatalError); } @@ -1048,14 +1036,14 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave setFaceInfo(changedFaces, changedFacesInfo); // Iterate until nothing changes - label iter = iterate(maxIter); + const label iter = iterate(maxIter); if ((maxIter > 0) && (iter >= maxIter)) { FatalErrorInFunction - << "Maximum number of iterations reached. Increase maxIter." << endl - << " maxIter:" << maxIter << endl - << " nChangedCells:" << changedCells_.size() << endl + << "Maximum number of iterations reached. Increase maxIter." << nl + << " maxIter:" << maxIter << nl + << " nChangedCells:" << changedCells_.size() << nl << " nChangedFaces:" << changedFaces_.size() << endl << exit(FatalError); } @@ -1065,14 +1053,14 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template<class Type, class TrackingData> -Foam::label Foam::FaceCellWave<Type, TrackingData>::getUnsetCells() const +Foam::label Foam::FaceCellWave<Type, TrackingData>::nUnvisitedCells() const { return nUnvisitedCells_; } template<class Type, class TrackingData> -Foam::label Foam::FaceCellWave<Type, TrackingData>::getUnsetFaces() const +Foam::label Foam::FaceCellWave<Type, TrackingData>::nUnvisitedFaces() const { return nUnvisitedFaces_; } @@ -1085,11 +1073,10 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::faceToCell() const labelList& owner = mesh_.faceOwner(); const labelList& neighbour = mesh_.faceNeighbour(); - label nInternalFaces = mesh_.nInternalFaces(); + const label nInternalFaces = mesh_.nInternalFaces(); - forAll(changedFaces_, changedFacei) + for (const label facei : changedFaces_) { - label facei = changedFaces_[changedFacei]; if (!changedFace_.test(facei)) { FatalErrorInFunction @@ -1098,42 +1085,43 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::faceToCell() << abort(FatalError); } - - const Type& neighbourWallInfo = allFaceInfo_[facei]; + const Type& newInfo = allFaceInfo_[facei]; // Evaluate all connected cells // Owner - label celli = owner[facei]; - Type& currentWallInfo = allCellInfo_[celli]; - - if (!currentWallInfo.equal(neighbourWallInfo, td_)) { - updateCell - ( - celli, - facei, - neighbourWallInfo, - propagationTol_, - currentWallInfo - ); + const label celli = owner[facei]; + Type& currInfo = allCellInfo_[celli]; + + if (!currInfo.equal(newInfo, td_)) + { + updateCell + ( + celli, + facei, + newInfo, + propagationTol_, + currInfo + ); + } } // Neighbour. if (facei < nInternalFaces) { - celli = neighbour[facei]; - Type& currentWallInfo2 = allCellInfo_[celli]; + const label celli = neighbour[facei]; + Type& currInfo = allCellInfo_[celli]; - if (!currentWallInfo2.equal(neighbourWallInfo, td_)) + if (!currInfo.equal(newInfo, td_)) { updateCell ( celli, facei, - neighbourWallInfo, + newInfo, propagationTol_, - currentWallInfo2 + currInfo ); } } @@ -1150,12 +1138,8 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::faceToCell() Pout<< " Changed cells : " << changedCells_.size() << endl; } - // Sum changedCells over all procs - label totNChanged = changedCells_.size(); - - reduce(totNChanged, sumOp<label>()); - - return totNChanged; + // Number of changedCells over all procs + return returnReduce(changedCells_.size(), sumOp<label>()); } @@ -1166,35 +1150,33 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::cellToFace() const cellList& cells = mesh_.cells(); - forAll(changedCells_, changedCelli) + for (const label celli : changedCells_) { - label celli = changedCells_[changedCelli]; - if (!changedCell_[celli]) + if (!changedCell_.test(celli)) { FatalErrorInFunction << "Cell " << celli << " not marked as having been changed" << abort(FatalError); } - const Type& neighbourWallInfo = allCellInfo_[celli]; + const Type& newInfo = allCellInfo_[celli]; // Evaluate all connected faces const labelList& faceLabels = cells[celli]; - forAll(faceLabels, faceLabelI) + for (const label facei : faceLabels) { - label facei = faceLabels[faceLabelI]; - Type& currentWallInfo = allFaceInfo_[facei]; + Type& currInfo = allFaceInfo_[facei]; - if (!currentWallInfo.equal(neighbourWallInfo, td_)) + if (!currInfo.equal(newInfo, td_)) { updateFace ( facei, celli, - neighbourWallInfo, + newInfo, propagationTol_, - currentWallInfo + currInfo ); } } @@ -1212,7 +1194,6 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::cellToFace() if (hasCyclicPatches_) { - // Transfer changed faces across cyclic halves handleCyclicPatches(); } @@ -1223,7 +1204,6 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::cellToFace() if (Pstream::parRun()) { - // Transfer changed faces from neighbouring processors. handleProcPatches(); } @@ -1232,12 +1212,9 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::cellToFace() Pout<< " Changed faces : " << changedFaces_.size() << endl; } - // Sum nChangedFaces over all procs - label totNChanged = changedFaces_.size(); - reduce(totNChanged, sumOp<label>()); - - return totNChanged; + // Number of changedFaces over all procs + return returnReduce(changedFaces_.size(), sumOp<label>()); } @@ -1245,9 +1222,13 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::cellToFace() template<class Type, class TrackingData> Foam::label Foam::FaceCellWave<Type, TrackingData>::iterate(const label maxIter) { + if (maxIter < 0) + { + return 0; + } + if (hasCyclicPatches_) { - // Transfer changed faces across cyclic halves handleCyclicPatches(); } @@ -1258,13 +1239,12 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::iterate(const label maxIter) if (Pstream::parRun()) { - // Transfer changed faces from neighbouring processors. handleProcPatches(); } label iter = 0; - while (iter < maxIter) + for (/*nil*/; iter < maxIter; ++iter) { if (debug) { @@ -1272,35 +1252,23 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::iterate(const label maxIter) } nEvals_ = 0; - - label nCells = faceToCell(); + const label nCells = faceToCell(); + const label nFaces = nCells ? cellToFace() : 0; if (debug) { - Info<< " Total changed cells : " << nCells << endl; + Info<< " Total evaluations : " + << nEvals_ << nl + << " Changed cells / faces : " + << nCells << " / " << nFaces << nl + << " Pending cells / faces : " + << nUnvisitedCells_ << " / " << nUnvisitedFaces_ << nl; } - if (nCells == 0) + if (!nCells || !nFaces) { break; } - - label nFaces = cellToFace(); - - if (debug) - { - Info<< " Total changed faces : " << nFaces << nl - << " Total evaluations : " << nEvals_ << nl - << " Remaining unvisited cells: " << nUnvisitedCells_ << nl - << " Remaining unvisited faces: " << nUnvisitedFaces_ << endl; - } - - if (nFaces == 0) - { - break; - } - - ++iter; } return iter; diff --git a/src/meshTools/algorithms/MeshWave/FaceCellWave.H b/src/meshTools/algorithms/MeshWave/FaceCellWave.H index e803a76aebf646efe61507b58077723c20d97817..29959f506944567ce29edfd3b79ebd87677bbcab 100644 --- a/src/meshTools/algorithms/MeshWave/FaceCellWave.H +++ b/src/meshTools/algorithms/MeshWave/FaceCellWave.H @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -57,7 +57,7 @@ SourceFiles namespace Foam { -// Forward declaration of classes +// Forward declarations class polyMesh; class polyPatch; @@ -79,16 +79,21 @@ class FaceCellWave { // Private Member Functions - //- Disallow default bitwise copy construct - FaceCellWave(const FaceCellWave&); + //- No copy construct + FaceCellWave(const FaceCellWave&) = delete; - //- Disallow default bitwise assignment - void operator=(const FaceCellWave&); + //- No copy assignment + void operator=(const FaceCellWave&) = delete; protected: - // Protected data + //- Information tagged with a source or destination id. + // With std::pair as lightweight, moveable container. + typedef std::pair<label,Type> taggedInfoType; + + + // Protected Data //- Reference to mesh const polyMesh& mesh_; @@ -108,15 +113,18 @@ protected: //- Has face changed bitSet changedFace_; - //- List of changed faces - DynamicList<label> changedFaces_; - //- Has cell changed bitSet changedCell_; + //- List of changed faces + DynamicList<label> changedFaces_; + // Cells that have changed DynamicList<label> changedCells_; + // Information exchange for explicit baffle connections + // Max capacity = 2x number of explicit connections + DynamicList<taggedInfoType> changedBaffles_; //- Contains cyclics const bool hasCyclicPatches_; @@ -132,8 +140,10 @@ protected: label nUnvisitedFaces_; - //- Updates cellInfo with information from neighbour. Updates all - // statistics. + // Protected Member Functions + + //- Updates cellInfo with information from neighbour. + // Updates all statistics. bool updateCell ( const label celli, @@ -143,8 +153,8 @@ protected: Type& cellInfo ); - //- Updates faceInfo with information from neighbour. Updates all - // statistics. + //- Updates faceInfo with information from neighbour. + // Updates all statistics. bool updateFace ( const label facei, @@ -154,8 +164,8 @@ protected: Type& faceInfo ); - //- Updates faceInfo with information from same face. Updates all - // statistics. + //- Updates faceInfo with information from same face. + // Updates all statistics. bool updateFace ( const label facei, @@ -179,8 +189,8 @@ protected: ( const polyPatch& patch, const label nFaces, - const labelList&, - const List<Type>& + const labelUList& changedFaces, + const List<Type>& changedFacesInfo ); //- Extract info for single patch only @@ -198,7 +208,7 @@ protected: ( const polyPatch& patch, const label nFaces, - const labelList& faceLabels, + const labelUList& faceLabels, List<Type>& faceInfo ) const; @@ -207,7 +217,7 @@ protected: ( const polyPatch& patch, const label nFaces, - const labelList& faceLabels, + const labelUList& faceLabels, List<Type>& faceInfo ) const; @@ -229,16 +239,18 @@ protected: ); //- Merge data from across processor boundaries + // Transfer changed faces from neighbouring processors. void handleProcPatches(); //- Merge data from across cyclics + // Transfer changed faces across cyclic halves void handleCyclicPatches(); //- Merge data from across AMI cyclics void handleAMICyclicPatches(); - //- Merge data across explicitly provided local connections (usually - // baffles) + //- Merge data across explicitly provided local connections + // These are usually baffles void handleExplicitConnections(); @@ -271,11 +283,11 @@ public: // Constructors - // Construct from mesh. Use setFaceInfo and iterate() to do actual - // calculation. + //- Construct from mesh. + //- Use setFaceInfo and iterate() to do actual calculation. FaceCellWave ( - const polyMesh&, + const polyMesh& mesh, UList<Type>& allFaceInfo, UList<Type>& allCellInfo, TrackingData& td = dummyTrackData_ @@ -283,11 +295,11 @@ public: //- Construct from mesh and list of changed faces with the Type // for these faces. Iterates until nothing changes or maxIter reached. - // (maxIter can be 0) + // (maxIter can be 0 or negative). 0 initializes, -1 does not FaceCellWave ( - const polyMesh&, - const labelList& initialChangedFaces, + const polyMesh& mesh, + const labelUList& initialChangedFaces, const List<Type>& changedFacesInfo, UList<Type>& allFaceInfo, UList<Type>& allCellInfo, @@ -298,13 +310,13 @@ public: //- Construct from mesh and explicitly connected boundary faces // and list of changed faces with the Type // for these faces. Iterates until nothing changes or maxIter reached. - // (maxIter can be 0) + // (maxIter can be 0 or negative). 0 initializes, -1 does not FaceCellWave ( - const polyMesh&, + const polyMesh& mesh, const labelPairList& explicitConnections, const bool handleCyclicAMI, - const labelList& initialChangedFaces, + const labelUList& initialChangedFaces, const List<Type>& changedFacesInfo, UList<Type>& allFaceInfo, UList<Type>& allCellInfo, @@ -314,70 +326,75 @@ public: //- Destructor - virtual ~FaceCellWave() - {} + virtual ~FaceCellWave() = default; // Member Functions - // Access + // Access + + //- Access allFaceInfo + UList<Type>& allFaceInfo() + { + return allFaceInfo_; + } + + //- Access allCellInfo + UList<Type>& allCellInfo() + { + return allCellInfo_; + } - //- Access allFaceInfo - UList<Type>& allFaceInfo() - { - return allFaceInfo_; - } + //- Additional data to be passed into container + const TrackingData& data() const + { + return td_; + } - //- Access allCellInfo - UList<Type>& allCellInfo() - { - return allCellInfo_; - } + //- Access mesh + const polyMesh& mesh() const + { + return mesh_; + } - //- Additional data to be passed into container - const TrackingData& data() const - { - return td_; - } + //- Get number of unvisited cells, + //- i.e. cells that were not (yet) reached from walking across mesh. + // This can happen from + // - not enough iterations done + // - a disconnected mesh + // - a mesh without walls in it + label nUnvisitedCells() const; - //- Access mesh - const polyMesh& mesh() const - { - return mesh_; - } + //- Get number of unvisited faces + label nUnvisitedFaces() const; - //- Get number of unvisited cells, i.e. cells that were not (yet) - // reached from walking across mesh. This can happen from - // - not enough iterations done - // - a disconnected mesh - // - a mesh without walls in it - label getUnsetCells() const; - //- Get number of unvisited faces - label getUnsetFaces() const; + // Edit + //- Set single initial changed face. + // This is a noop if the face had already been visited + void setFaceInfo(const label facei, const Type& faceInfo); - // Edit + //- Set initial changed faces + void setFaceInfo + ( + const labelUList& changedFaces, + const List<Type>& changedFacesInfo + ); - //- Set initial changed faces - void setFaceInfo - ( - const labelList& changedFaces, - const List<Type>& changedFacesInfo - ); + //- Propagate from face to cell. + // \return total number of cells (over all processors) changed. + virtual label faceToCell(); - //- Propagate from face to cell. Returns total number of cells - // (over all processors) changed. - virtual label faceToCell(); + //- Propagate from cell to face. + // \return total number of faces (over all processors) changed. + // Note that faces on processor patches are counted twice. + virtual label cellToFace(); - //- Propagate from cell to face. Returns total number of faces - // (over all processors) changed. (Faces on processorpatches are - // counted double) - virtual label cellToFace(); + //- Iterate until no changes or maxIter reached. + // \return the number of iterations taken. + virtual label iterate(const label maxIter); - //- Iterate until no changes or maxIter reached. Returns actual - // number of iterations. - virtual label iterate(const label maxIter); }; diff --git a/src/meshTools/algorithms/MeshWave/MeshWave.H b/src/meshTools/algorithms/MeshWave/MeshWave.H index a45ba4b9884e1fc4ef5c02eb157d1dfc5123dd60..f522a30bd05b43e35332bd8fda566814ced99baa 100644 --- a/src/meshTools/algorithms/MeshWave/MeshWave.H +++ b/src/meshTools/algorithms/MeshWave/MeshWave.H @@ -141,20 +141,20 @@ public: return calc_.iterate(maxIter); } - //- Get number of unvisited cells, i.e. cells that were not (yet) + //- Number of unvisited cells, i.e. cells that were not (yet) // reached from walking across mesh. This can happen from // - not enough iterations done // - a disconnected mesh // - a mesh without walls in it - label getUnsetCells() const + label nUnvisitedCells() const { - return calc_.getUnsetCells(); + return calc_.nUnvisitedCells(); } - //- Get number of unvisited faces - label getUnsetFaces() const + //- Number of unvisited faces + label nUnvisitedFaces() const { - return calc_.getUnsetFaces(); + return calc_.nUnvisitedFaces(); } }; diff --git a/src/meshTools/algorithms/PatchEdgeFaceWave/PatchEdgeFaceWave.C b/src/meshTools/algorithms/PatchEdgeFaceWave/PatchEdgeFaceWave.C index 24a30b6864c3c2653f49bbda5aae01672744d5a7..7e3376cb5be62e717dc54c434103d94e61c4b680 100644 --- a/src/meshTools/algorithms/PatchEdgeFaceWave/PatchEdgeFaceWave.C +++ b/src/meshTools/algorithms/PatchEdgeFaceWave/PatchEdgeFaceWave.C @@ -420,7 +420,7 @@ template class TrackingData > Foam::label Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>:: -getUnsetEdges() const +nUnvisitedEdges() const { return nUnvisitedEdges_; } @@ -433,7 +433,7 @@ template class TrackingData > Foam::label Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>:: -getUnsetFaces() const +nUnvisitedFaces() const { return nUnvisitedFaces_; } diff --git a/src/meshTools/algorithms/PatchEdgeFaceWave/PatchEdgeFaceWave.H b/src/meshTools/algorithms/PatchEdgeFaceWave/PatchEdgeFaceWave.H index 67ee9362f7113fd095de2505d03bb9accfdb2a25..176aaded9d7a21313937c27dadbd9410cafd7ee8 100644 --- a/src/meshTools/algorithms/PatchEdgeFaceWave/PatchEdgeFaceWave.H +++ b/src/meshTools/algorithms/PatchEdgeFaceWave/PatchEdgeFaceWave.H @@ -226,14 +226,14 @@ public: return td_; } - //- Get number of unvisited faces, i.e. faces that were not (yet) + //- Number of unvisited faces, i.e. faces that were not (yet) // reached from walking across patch. This can happen from // - not enough iterations done // - a disconnected patch // - a patch without walls in it - label getUnsetFaces() const; + label nUnvisitedFaces() const; - label getUnsetEdges() const; + label nUnvisitedEdges() const; //- Copy initial data into allEdgeInfo_ void setEdgeInfo diff --git a/src/meshTools/algorithms/PointEdgeWave/PointEdgeWave.C b/src/meshTools/algorithms/PointEdgeWave/PointEdgeWave.C index 668f1f9b6710eae1bd13235c917d646c425df9c9..435b59d2f18963b93f4828d718fca931ffc76c1b 100644 --- a/src/meshTools/algorithms/PointEdgeWave/PointEdgeWave.C +++ b/src/meshTools/algorithms/PointEdgeWave/PointEdgeWave.C @@ -714,14 +714,14 @@ Foam::PointEdgeWave<Type, TrackingData>::~PointEdgeWave() template<class Type, class TrackingData> -Foam::label Foam::PointEdgeWave<Type, TrackingData>::getUnsetPoints() const +Foam::label Foam::PointEdgeWave<Type, TrackingData>::nUnvisitedPoints() const { return nUnvisitedPoints_; } template<class Type, class TrackingData> -Foam::label Foam::PointEdgeWave<Type, TrackingData>::getUnsetEdges() const +Foam::label Foam::PointEdgeWave<Type, TrackingData>::nUnvisitedEdges() const { return nUnvisitedEdges_; } @@ -737,9 +737,9 @@ void Foam::PointEdgeWave<Type, TrackingData>::setPointInfo { forAll(changedPoints, changedPointi) { - label pointi = changedPoints[changedPointi]; + const label pointi = changedPoints[changedPointi]; - bool wasValid = allPointInfo_[pointi].valid(td_); + const bool wasValid = allPointInfo_[pointi].valid(td_); // Copy info for pointi allPointInfo_[pointi] = changedPointsInfo[changedPointi]; diff --git a/src/meshTools/algorithms/PointEdgeWave/PointEdgeWave.H b/src/meshTools/algorithms/PointEdgeWave/PointEdgeWave.H index 6da80f46e40b65359283947037d868e9e79668c9..ad85976f01271e828333189922e097600f0e7ccf 100644 --- a/src/meshTools/algorithms/PointEdgeWave/PointEdgeWave.H +++ b/src/meshTools/algorithms/PointEdgeWave/PointEdgeWave.H @@ -287,14 +287,14 @@ public: return td_; } - //- Get number of unvisited edges, i.e. edges that were not (yet) + //- Number of unvisited edges, i.e. edges that were not (yet) // reached from walking across mesh. This can happen from // - not enough iterations done // - a disconnected mesh // - a mesh without walls in it - label getUnsetEdges() const; + label nUnvisitedEdges() const; - label getUnsetPoints() const; + label nUnvisitedPoints() const; //- Copy initial data into allPointInfo_ void setPointInfo diff --git a/src/meshTools/regionSplit/regionSplit.C b/src/meshTools/regionSplit/regionSplit.C index a901e8d763b054ddbb7c3dd3bfd6c50929989e81..eb399eae966692dbd942eb168ce8cfe69ac161c1 100644 --- a/src/meshTools/regionSplit/regionSplit.C +++ b/src/meshTools/regionSplit/regionSplit.C @@ -2,8 +2,8 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation - \\/ M anipulation | Copyright (C) 2015-2018 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\/ M anipulation | Copyright (C) 2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -24,12 +24,11 @@ License \*---------------------------------------------------------------------------*/ #include "regionSplit.H" -#include "FaceCellWave.H" #include "cyclicPolyPatch.H" #include "processorPolyPatch.H" #include "globalIndex.H" #include "syncTools.H" -#include "minData.H" +#include "clockValue.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -38,443 +37,396 @@ namespace Foam defineTypeNameAndDebug(regionSplit, 0); } -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // +static constexpr Foam::label UNASSIGNED = -1; +static constexpr Foam::label BLOCKED = -2; -void Foam::regionSplit::calcNonCompactRegionSplit -( - const globalIndex& globalFaces, - const boolList& blockedFace, - const List<labelPair>& explicitConnections, - labelList& cellRegion -) const +// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // + +namespace Foam { - // Field on cells and faces. - List<minData> cellData(mesh().nCells()); - List<minData> faceData(mesh().nFaces()); - // Take over blockedFaces by seeding a negative number - // (so is always less than the decomposition) +//- Use key/value from mapper to modify the values for map + +template<class Container> +static void inplaceMapValue(const Map<label>& mapper, Container& input) +{ + if (mapper.empty()) + { + return; + } - forAll(blockedFace, facei) + for (auto iter = input.begin(); iter != input.end(); ++iter) { - if (blockedFace[facei]) + label& value = *iter; + + auto mapIter = mapper.find(value); + if (mapIter.found()) { - faceData[facei] = minData(-2); + value = *mapIter; } } +} - // Seed all faces on (real) boundaries and cell faces next to blockFace, - // since regions can only occur because of boundaries (or blocked faces) - bitSet isSeed(mesh().nFaces()); +//- The sizes of a List of containers (eg, labelHashSet) +template<class Container> +static labelList containerSizes(const UList<Container>& input) +{ + const label len = input.size(); - // Get internal or coupled faces - bitSet isConnection(syncTools::getInternalOrCoupledFaces(mesh())); + labelList output(len); - // 1. Seed (real) boundaries - for - ( - label facei = mesh().nInternalFaces(); - facei < mesh().nFaces(); - ++facei - ) + for (label i = 0; i < len; ++i) { - if (!isConnection[facei]) - { - isSeed.set(facei); - } + output[i] = input[i].size(); } - // 2. Seed (internal) faces on cells next to blocked (internal) faces - // (since we can't seed the blocked faces themselves) + return output; +} - if (blockedFace.size()) - { - label nBlockedCells = 0; +} // end of namespace Foam - label celli = 0; - for (const cell& cFaces : mesh().cells()) - { - bool blockedCell = false; - label connectedFacei = -1; - for (const label facei : cFaces) - { - if (blockedFace[facei]) - { - blockedCell = true; - } - else if (isConnection[facei]) - { - connectedFacei = facei; - } - } - - if (blockedCell) - { - if (connectedFacei < 0) - { - ++nBlockedCells; +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - if (debug && nBlockedCells <= 10) - { - // Report a few, but not all - WarningInFunction - << "Cell " << celli << " has all faces blocked." - << endl; - } - } - else - { - isSeed.set(connectedFacei); - } - } +void Foam::regionSplit::updateFacePair +( + const label face0, + const label face1, - ++celli; + labelList& faceRegion, + DynamicList<label>& facesChanged +) const +{ + if (faceRegion[face0] == UNASSIGNED) + { + if (faceRegion[face1] >= 0) + { + faceRegion[face0] = faceRegion[face1]; + facesChanged.append(face0); } - - reduce(nBlockedCells, sumOp<label>()); - - if (nBlockedCells) + else if (faceRegion[face1] == BLOCKED) { - Info<< "Detected " << nBlockedCells << " from " - << returnReduce(mesh().nCells(), sumOp<label>()) - << " cells with all faces blocked." << nl; + // face1 blocked but not face0. + // - illegal for coupled faces, OK for explicit connections. } } - - List<label> seedFaces(isSeed.toc()); - List<minData> seedData(seedFaces.size()); - - // Seed face with globally unique number - forAll(seedFaces, i) + else if (faceRegion[face0] >= 0) { - const label facei = seedFaces[i]; - seedData[i] = minData(globalFaces.toGlobal(facei)); - } - - if (debug) - { - Pout<< "Seeded " << seedFaces.size() - << " boundary and internal faces out of" - << " total:" << mesh().nInternalFaces() - << " boundary:" << mesh().nFaces()-mesh().nInternalFaces() - << endl; - } - - // Propagate information inwards - FaceCellWave<minData> deltaCalc - ( - mesh(), - explicitConnections, - false, // disable walking through cyclicAMI for backwards compatibility - seedFaces, - seedData, - faceData, - cellData, - mesh().globalData().nTotalCells()+1 - ); - - - // And extract - cellRegion.setSize(mesh().nCells()); - forAll(cellRegion, celli) - { - if (cellData[celli].valid(deltaCalc.data())) + if (faceRegion[face1] == UNASSIGNED) { - cellRegion[celli] = cellData[celli].data(); + faceRegion[face1] = faceRegion[face0]; + facesChanged.append(face1); } - else + else if (faceRegion[face1] == BLOCKED) { - // Unvisited cell -> only possible if surrounded by blocked faces - // since: - // - all walls become seed faces - // - all faces on cells with blocked faces become seed faces - // If so make up region from any of the faces - const cell& cFaces = mesh().cells()[celli]; - label facei = cFaces[0]; - - if (blockedFace.size() && !blockedFace[facei]) - { - FatalErrorInFunction - << "Problem: unblocked face " << facei - << " at " << mesh().faceCentres()[facei] - << " on unassigned cell " << celli - << mesh().cellCentres()[celli] - << exit(FatalError); - } - cellRegion[celli] = globalFaces.toGlobal(facei); + // face1 blocked but not face0. + // - illegal for coupled faces, OK for explicit connections. + } + else if (faceRegion[face1] != faceRegion[face0]) + { + FatalErrorInFunction + << "Problem : coupled face " << face0 + << " on patch " << mesh().boundaryMesh().whichPatch(face0) + << " has region " << faceRegion[face0] + << " but coupled face " << face1 + << " has region " << faceRegion[face1] << nl + << "Is your blocked faces specification" + << " synchronized across coupled boundaries?" << endl + << abort(FatalError); } } } -Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit +void Foam::regionSplit::fillSeedMask ( - const bool doGlobalRegions, - const boolList& blockedFace, const List<labelPair>& explicitConnections, - - labelList& cellRegion + const label seedCellId, + const label markValue, + labelList& cellRegion, + labelList& faceRegion ) const { - // See header in regionSplit.H + // Seed cell + cellRegion[seedCellId] = markValue; + // Faces on seed cell + changedFaces_.clear(); + for (const label facei : mesh().cells()[seedCellId]) + { + if (faceRegion[facei] == UNASSIGNED) + { + faceRegion[facei] = markValue; + changedFaces_.append(facei); + } + } - if (!doGlobalRegions) + const polyBoundaryMesh& patches = mesh().boundaryMesh(); + + // Loop over changed faces. FaceCellWave in small. + + while (changedFaces_.size()) { - // Block all parallel faces to avoid comms across - boolList coupledOrBlockedFace(blockedFace); - const polyBoundaryMesh& pbm = mesh().boundaryMesh(); + changedCells_.clear(); - if (coupledOrBlockedFace.size()) + for (const label facei : changedFaces_) { - forAll(pbm, patchi) + const label own = mesh().faceOwner()[facei]; + + if (cellRegion[own] == UNASSIGNED) + { + cellRegion[own] = markValue; + changedCells_.append(own); + } + + if (mesh().isInternalFace(facei)) { - const polyPatch& pp = pbm[patchi]; - if (isA<processorPolyPatch>(pp)) + const label nei = mesh().faceNeighbour()[facei]; + + if (cellRegion[nei] == UNASSIGNED) { - label facei = pp.start(); - forAll(pp, i) - { - coupledOrBlockedFace[facei++] = true; - } + cellRegion[nei] = markValue; + changedCells_.append(nei); } } } - // Create dummy (local only) globalIndex - labelList offsets(Pstream::nProcs()+1, 0); - for (label i = Pstream::myProcNo()+1; i < offsets.size(); i++) + if (debug & 2) { - offsets[i] = mesh().nFaces(); + Pout<< " Changed cells / faces : " + << changedCells_.size() << " / " << changedFaces_.size() + << " before sync" << endl; } - const globalIndex globalRegions(std::move(offsets)); - - // Minimise regions across connected cells - // Note: still uses global decisions so all processors are running - // in lock-step, i.e. slowest determines overall time. - // To avoid this we could switch off Pstream::parRun. - calcNonCompactRegionSplit - ( - globalRegions, - coupledOrBlockedFace, - explicitConnections, - cellRegion - ); - // Compact - Map<label> globalToCompact(mesh().nCells()/8); - forAll(cellRegion, celli) + // Loop over changedCells and collect faces + changedFaces_.clear(); + for (const label celli : changedCells_) { - label region = cellRegion[celli]; - - label globalRegion; - - const auto fnd = globalToCompact.cfind(region); - if (fnd.found()) - { - globalRegion = *fnd; - } - else + for (const label facei : mesh().cells()[celli]) { - globalRegion = globalRegions.toGlobal(globalToCompact.size()); - globalToCompact.insert(region, globalRegion); + if (faceRegion[facei] == UNASSIGNED) + { + faceRegion[facei] = markValue; + changedFaces_.append(facei); + } } - cellRegion[celli] = globalRegion; } - // Return globalIndex with size = localSize and all regions local - labelList compactOffsets(Pstream::nProcs()+1, 0); - for (label i = Pstream::myProcNo()+1; i < compactOffsets.size(); i++) + // Update locally coupled faces + // Global connections are done later. + + for (const polyPatch& pp : patches) { - compactOffsets[i] = globalToCompact.size(); - } + if + ( + isA<cyclicPolyPatch>(pp) + && refCast<const cyclicPolyPatch>(pp).owner() + ) + { + // Transfer from neighbourPatch to here or vice versa. - return autoPtr<globalIndex>::New(std::move(compactOffsets)); - } + const cyclicPolyPatch& cycPatch = + refCast<const cyclicPolyPatch>(pp); + label face0 = cycPatch.start(); + forAll(cycPatch, i) + { + const label face1 = cycPatch.transformGlobalFace(face0); - // Initial global region numbers - const globalIndex globalRegions(mesh().nFaces()); + updateFacePair + ( + face0, + face1, + faceRegion, + changedFaces_ + ); - // Minimise regions across connected cells (including parallel) - calcNonCompactRegionSplit - ( - globalRegions, - blockedFace, - explicitConnections, - cellRegion - ); + ++face0; + } + } + } + for (const labelPair& pr : explicitConnections) + { + updateFacePair + ( + pr.first(), + pr.second(), + faceRegion, + changedFaces_ + ); + } - // Now our cellRegion will have - // - non-local regions (i.e. originating from other processors) - // - non-compact locally originating regions - // so we'll need to compact + if (debug & 2) + { + Pout<< " Changed faces : " + << changedFaces_.size() + << " after sync" << endl; + } + } +} - // 4a: count per originating processor the number of regions - labelList nOriginating(Pstream::nProcs(), 0); - { - labelHashSet haveRegion(mesh().nCells()/8); - forAll(cellRegion, celli) +Foam::label Foam::regionSplit::calcLocalRegionSplit +( + const boolList& blockedFace, + const List<labelPair>& explicitConnections, + + labelList& cellRegion +) const +{ + clockValue timing(debug); + + if (debug) + { + if (blockedFace.size()) { - label region = cellRegion[celli]; + // Check that blockedFace is synced. + boolList syncBlockedFace(blockedFace); + syncTools::swapFaceList(mesh(), syncBlockedFace); - // Count originating processor. Use isLocal as efficiency since - // most cells are locally originating. - if (globalRegions.isLocal(region)) - { - if (haveRegion.insert(region)) - { - nOriginating[Pstream::myProcNo()]++; - } - } - else + forAll(syncBlockedFace, facei) { - label proci = globalRegions.whichProcID(region); - if (haveRegion.insert(region)) + if (syncBlockedFace[facei] != blockedFace[facei]) { - nOriginating[proci]++; + FatalErrorInFunction + << "Face " << facei << " not synchronised. My value:" + << blockedFace[facei] << " coupled value:" + << syncBlockedFace[facei] + << abort(FatalError); } } } } - if (debug) - { - Pout<< "Counted " << nOriginating[Pstream::myProcNo()] - << " local regions." << endl; - } + changedCells_.reserve(mesh_.nCells()); + changedFaces_.reserve(mesh_.nFaces()); + // Region per face. + // -1 = unassigned + // -2 = blocked + labelList faceRegion(mesh().nFaces(), UNASSIGNED); - // Global numbering for compacted local regions - autoPtr<globalIndex> globalCompactPtr - ( - new globalIndex(nOriginating[Pstream::myProcNo()]) - ); - const globalIndex& globalCompact = globalCompactPtr(); - - - // 4b: renumber - // Renumber into compact indices. Note that since we've already made - // all regions global we now need a Map to store the compacting information - // instead of a labelList - otherwise we could have used a straight - // labelList. - - // Local compaction map - Map<label> globalToCompact(2*nOriginating[Pstream::myProcNo()]); - // Remote regions we want the compact number for - List<labelHashSet> nonLocal(Pstream::nProcs()); - forAll(nonLocal, proci) + if (blockedFace.size()) { - if (proci != Pstream::myProcNo()) + forAll(blockedFace, facei) { - nonLocal[proci].resize(2*nOriginating[proci]); + if (blockedFace[facei]) + { + faceRegion[facei] = BLOCKED; + } } } - forAll(cellRegion, celli) + + // Assign local regions + // ~~~~~~~~~~~~~~~~~~~~ + + // Start with region 0 + label nLocalRegions = 0; + + for (label seedCellId = 0; seedCellId < cellRegion.size(); ++seedCellId) { - label region = cellRegion[celli]; - if (globalRegions.isLocal(region)) + // Find next unset cell - use as seed + + for (; seedCellId < cellRegion.size(); ++seedCellId) { - // Insert new compact region (if not yet present) - globalToCompact.insert - ( - region, - globalCompact.toGlobal(globalToCompact.size()) - ); + if (cellRegion[seedCellId] == UNASSIGNED) + { + break; + } } - else + + if (seedCellId >= cellRegion.size()) { - nonLocal[globalRegions.whichProcID(region)].insert(region); + break; } - } - - // Now we have all the local regions compacted. Now we need to get the - // non-local ones from the processors to whom they are local. - // Convert the nonLocal (labelHashSets) to labelLists. + fillSeedMask + ( + explicitConnections, + seedCellId, + nLocalRegions, + cellRegion, + faceRegion + ); - labelListList sendNonLocal(Pstream::nProcs()); - forAll(sendNonLocal, proci) - { - sendNonLocal[proci] = nonLocal[proci].toc(); + ++nLocalRegions; // Next region } + // Discard temporary working data + changedCells_.clearStorage(); + changedFaces_.clearStorage(); + if (debug) { - forAll(sendNonLocal, proci) + forAll(cellRegion, celli) { - Pout<< " from processor " << proci - << " want " << sendNonLocal[proci].size() - << " region numbers." - << endl; + if (cellRegion[celli] < 0) + { + FatalErrorInFunction + << "cell:" << celli << " region:" << cellRegion[celli] + << abort(FatalError); + } } - Pout<< endl; - } - - - // Get the wanted region labels into recvNonLocal - labelListList recvNonLocal; - Pstream::exchange<labelList, label>(sendNonLocal, recvNonLocal); - - // Now we have the wanted compact region labels that proci wants in - // recvNonLocal[proci]. Construct corresponding list of compact - // region labels to send back. - labelListList sendWantedLocal(Pstream::nProcs()); - forAll(recvNonLocal, proci) - { - const labelList& nonLocal = recvNonLocal[proci]; - sendWantedLocal[proci].setSize(nonLocal.size()); - - forAll(nonLocal, i) + forAll(faceRegion, facei) { - sendWantedLocal[proci][i] = globalToCompact[nonLocal[i]]; + if (faceRegion[facei] == UNASSIGNED) + { + FatalErrorInFunction + << "face:" << facei << " region:" << faceRegion[facei] + << abort(FatalError); + } } } + if (debug) + { + Info<<"regionSplit = " << double(timing.elapsed()) << "s\n"; + } - // Send back (into recvNonLocal) - recvNonLocal.clear(); - Pstream::exchange<labelList, label>(sendWantedLocal, recvNonLocal); - sendWantedLocal.clear(); + return nLocalRegions; +} - // Now recvNonLocal contains for every element in setNonLocal the - // corresponding compact number. Insert these into the local compaction - // map. - forAll(recvNonLocal, proci) - { - const labelList& wantedRegions = sendNonLocal[proci]; - const labelList& compactRegions = recvNonLocal[proci]; +Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit +( + const bool doGlobalRegions, + const boolList& blockedFace, + const List<labelPair>& explicitConnections, - forAll(wantedRegions, i) - { - globalToCompact.insert(wantedRegions[i], compactRegions[i]); - } - } + labelList& cellRegion +) const +{ + const label nLocalRegions = calcLocalRegionSplit + ( + blockedFace, + explicitConnections, + cellRegion + ); - // Finally renumber the regions - forAll(cellRegion, celli) + if (!doGlobalRegions) { - cellRegion[celli] = globalToCompact[cellRegion[celli]]; + return autoPtr<globalIndex>::New(nLocalRegions); } - return globalCompactPtr; + return reduceRegions(nLocalRegions, blockedFace, cellRegion); } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::regionSplit::regionSplit(const polyMesh& mesh, const bool doGlobalRegions) +Foam::regionSplit::regionSplit +( + const polyMesh& mesh, + const bool doGlobalRegions +) : MeshObject<polyMesh, Foam::TopologicalMeshObject, regionSplit>(mesh), labelList(mesh.nCells(), -1) @@ -530,4 +482,354 @@ Foam::regionSplit::regionSplit } +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::reduceRegions +( + const label numLocalRegions, + const boolList& blockedFace, + + labelList& cellRegion +) const +{ + clockValue timing(debug); + + if (cellRegion.size() != mesh().nCells()) + { + FatalErrorInFunction + << "The cellRegion size " << cellRegion.size() + << " is not equal to the of number of cells " + << mesh().nCells() << endl + << abort(FatalError); + } + + + // (numLocalRegions < 0) to signal that region information should be + // determined ourselves. This is not really efficient, but can be useful + + const label nLocalRegions = + ( + numLocalRegions < 0 + ? labelHashSet(cellRegion).size() + : numLocalRegions + ); + + + // Preliminary global region numbers + const globalIndex globalRegions(nLocalRegions); + + + // Lookup table of local region to global region. + // Initially an identity mapping of the uncombined global values. + + Map<label> localToGlobal(2*nLocalRegions); + for (const label regioni : cellRegion) + { + localToGlobal.insert(regioni, globalRegions.toGlobal(regioni)); + } + + // To update the localToGlobal mapping during traversal of the boundaries + // and later when finalizing things. + Map<label> updateLookup(2*nLocalRegions); + + + // Note that we use two separate maps during the process. + // The localToGlobal is used to map the local to global regions. + // Merging across processors will normally make this a many->few mapping. + // However, we may need to walk up and down processor boundaries several + // times before all the information propagates through. + // During these traversals, it will normally be more efficient to just + // update the mapping without updating the cellRegion immediately. + // Only after everything is finalized do we renumber all of the cell + // regions. + + + // Merge global regions + // ~~~~~~~~~~~~~~~~~~~~ + // Regions across non-blocked proc patches get merged. + // This will set merged global regions to be the min of both. + // (this will create gaps in the global region list so they will get + // merged later on) + + const polyBoundaryMesh& patches = mesh().boundaryMesh(); + + // Buffer for swapping boundary information + labelList nbrRegion(mesh().nFaces()-mesh().nInternalFaces()); + + do + { + if (debug) + { + Pout<< nl << "-- Starting Iteration --" << endl; + } + + updateLookup.clear(); + nbrRegion = UNASSIGNED; + + // Region information to send + for (const polyPatch& pp : patches) + { + if (pp.coupled()) + { + const labelUList& faceCells = pp.faceCells(); + SubList<label> patchNbrRegion + ( + nbrRegion, + pp.size(), + pp.start()-mesh().nInternalFaces() + ); + + forAll(faceCells, patchFacei) + { + const label meshFacei = pp.start()+patchFacei; + + if (!blockedFace[meshFacei]) + { + // Send the most currently updated region Id + const label orig = cellRegion[faceCells[patchFacei]]; + + patchNbrRegion[patchFacei] = localToGlobal[orig]; + } + } + } + } + syncTools::swapBoundaryFaceList(mesh(), nbrRegion); + + // Receive and reduce region information + for (const polyPatch& pp : patches) + { + if (pp.coupled()) + { + const labelUList& faceCells = pp.faceCells(); + SubList<label> patchNbrRegion + ( + nbrRegion, + pp.size(), + pp.start()-mesh().nInternalFaces() + ); + + forAll(faceCells, patchFacei) + { + const label meshFacei = pp.start()+patchFacei; + + if (!blockedFace[meshFacei]) + { + // Reduction by retaining the min region id. + + const label orig = cellRegion[faceCells[patchFacei]]; + + const label sent = localToGlobal[orig]; + const label recv = patchNbrRegion[patchFacei]; + + // Record the minimum value seen + if (recv < sent) + { + auto fnd = updateLookup.find(sent); + if (!fnd.found()) + { + updateLookup.insert(sent, recv); + } + else if (recv < *fnd) + { + *fnd = recv; + } + } + } + } + } + } + + + // Note: by always using the minimum region number across the + // processor faces, we effect a consolidation of connected regions + // and converge to a unique number for each distinct region. + + + // Update localToGlobal according to newly exchanged information + + inplaceMapValue(updateLookup, localToGlobal); + + if (debug & 2) + { + labelList keys(localToGlobal.sortedToc()); + labelList vals(keys.size()); + forAll(keys, i) + { + vals[i] = localToGlobal[keys[i]]; + } + + Pout<< "Updated local regions:" << nl + << "old: " << flatOutput(keys) << nl + << "new: " << flatOutput(vals) << endl; + } + else if (debug) + { + Pout<< "Updated " << localToGlobal.size() + << " local regions" << endl; + } + + // Continue until there are no further changes + } + while (returnReduce(!updateLookup.empty(), orOp<bool>())); + + + // + // We will need compact numbers locally and non-locally + // + + // Determine the local compact numbering + label nCompact = 0; + { + labelHashSet localRegion(2*localToGlobal.size()); + + forAllConstIters(localToGlobal, iter) + { + const label regioni = iter.object(); + + if (globalRegions.isLocal(regioni)) + { + localRegion.insert(regioni); + } + } + + nCompact = localRegion.size(); + } + + + // The new global numbering using compacted local regions + auto globalCompactPtr = autoPtr<globalIndex>::New(nCompact); + const auto& globalCompact = *globalCompactPtr; + + + // Determine the following: + // - the local compact regions (store as updateLookup) + // - the non-local regions, ordered according to the processor on which + // they are local. + + + // The local compaction map (updated local to compact local numbering) + updateLookup.clear(); + + labelListList sendNonLocal(Pstream::nProcs()); + + { + List<labelHashSet> nonLocal(Pstream::nProcs(), labelHashSet(0)); + + // Use estimate of sizing for non-local regions + forAll(nonLocal, proci) + { + if (proci != Pstream::myProcNo()) + { + nonLocal[proci].resize + ( + 2*((nLocalRegions-nCompact)/Pstream::nProcs()) + ); + } + } + + + forAllConstIters(localToGlobal, iter) + { + const label regioni = iter.object(); + + if (globalRegions.isLocal(regioni)) + { + updateLookup.insert + ( + regioni, + globalCompact.toGlobal(updateLookup.size()) + ); + } + else + { + nonLocal[globalRegions.whichProcID(regioni)].insert(regioni); + } + } + + if (debug) + { + Pout<< " per processor nonLocal regions: " + << flatOutput(containerSizes(nonLocal)) << endl; + } + + + // Convert to label list + forAll(sendNonLocal, proci) + { + sendNonLocal[proci] = nonLocal[proci].toc(); + } + } + + + // Get the wanted region labels into recvNonLocal + labelListList recvNonLocal; + Pstream::exchange<labelList, label> + ( + sendNonLocal, + recvNonLocal + ); + + + // The recvNonLocal[proci] region labels are what proci requires. + // Transcribe into their compacted number. + + { + labelListList sendLocal(std::move(recvNonLocal)); + + for (labelList& send : sendLocal) + { + for (label& regioni : send) + { + regioni = updateLookup[regioni]; + } + } + + // Send back (into recvNonLocal) + Pstream::exchange<labelList, label> + ( + sendLocal, + recvNonLocal + ); + } + + + // Now recvNonLocal and sendNonLocal contain matched pairs with + // sendNonLocal being the non-compact region and recvNonLocal being + // the compact region. + // + // Insert these into the local compaction map. + + forAll(recvNonLocal, proci) + { + const labelList& send = sendNonLocal[proci]; + const labelList& recv = recvNonLocal[proci]; + + forAll(send, i) + { + updateLookup.insert(send[i], recv[i]); + } + } + + + // Now renumber the localToGlobal to use the final compact global values + inplaceMapValue(updateLookup, localToGlobal); + + + // Can now finally use localToGlobal to renumber cellRegion + + forAll(cellRegion, celli) + { + cellRegion[celli] = localToGlobal[cellRegion[celli]]; + } + + if (debug) + { + Info<<"regionSplit::reduceRegions = " + << double(timing.elapsed()) << "s\n"; + } + + return globalCompactPtr; +} + + // ************************************************************************* // diff --git a/src/meshTools/regionSplit/regionSplit.H b/src/meshTools/regionSplit/regionSplit.H index 1051491f190698576abd531981d80a2e49768b0c..494e69cbf50b92ab546395a4784a548862e509dd 100644 --- a/src/meshTools/regionSplit/regionSplit.H +++ b/src/meshTools/regionSplit/regionSplit.H @@ -2,8 +2,8 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation - \\/ M anipulation | + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\/ M anipulation | Copyright (C) 2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -103,12 +103,11 @@ Description \endverbatim - Can optionally keep all regions local to the processor. - Note: does not walk across cyclicAMI/cyclicACMI - since not 'coupled()' - at the patch level. - +Note + does not walk across cyclicAMI/cyclicACMI - since these are not + \c coupled() at the patch level. SourceFiles regionSplit.C @@ -128,6 +127,8 @@ SourceFiles namespace Foam { +// Forward declarations + class polyMesh; /*---------------------------------------------------------------------------*\ @@ -141,21 +142,52 @@ class regionSplit { // Private data - mutable autoPtr<globalIndex> globalNumberingPtr_; + autoPtr<globalIndex> globalNumberingPtr_; + + //- Temporary list of cells that have changed + mutable DynamicList<label> changedCells_; + + //- Temporary list of faces that have changed + mutable DynamicList<label> changedFaces_; // Private Member Functions - //- Calculate region split in non-compact (global) numbering. - void calcNonCompactRegionSplit + //- Update faceRegion data between (non-processor) coupled faces. + void updateFacePair + ( + const label face0, + const label face1, + labelList& faceRegion, + DynamicList<label>& facesChanged + ) const; + + //- Given a seed cell label, fill cellRegion/faceRegion with markValue + //- for contiguous region around it + void fillSeedMask + ( + const List<labelPair>& explicitConnections, + const label seedCellID, + const label markValue, + labelList& cellRegion, + labelList& faceRegion + ) const; + + + //- Calculate the local region split. + // \return number of processor-local regions, + // without consolidation between procesors + label calcLocalRegionSplit ( - const globalIndex& globalFaces, const boolList& blockedFace, const List<labelPair>& explicitConnections, labelList& cellRegion ) const; - //- Calculate global region split. Return globalIndex. + + //- Calculate the local region split. + // \return number of processor-local regions, + // without consolidation between procesors autoPtr<globalIndex> calcRegionSplit ( const bool doGlobalRegions, @@ -176,27 +208,27 @@ public: //- Construct from mesh regionSplit ( - const polyMesh&, + const polyMesh& mesh, const bool doGlobalRegions = Pstream::parRun() ); //- Construct from mesh and whether face is blocked - // \note blockedFace has to be consistent across coupled faces! + // \note blockedFace must be consistent across coupled faces! regionSplit ( - const polyMesh&, + const polyMesh& mesh, const boolList& blockedFace, const bool doGlobalRegions = Pstream::parRun() ); - //- Construct from mesh and whether face is blocked. - // Additional explicit connections between normal boundary faces. - // \note blockedFace has to be consistent across coupled faces! + //- Construct from mesh and whether face is blocked, with additional + //- explicit connections between normal boundary faces. + // \note blockedFace must be consistent across coupled faces! regionSplit ( - const polyMesh&, + const polyMesh& mesh, const boolList& blockedFace, - const List<labelPair>&, + const List<labelPair>& explicitConnections, const bool doGlobalRegions = Pstream::parRun() ); @@ -220,6 +252,19 @@ public: { return globalNumbering().size(); } + + + //- Manually consolidate the regions globally by swapping information + // between processor domains and reducing the regions accordingly. + // + // \return number of local regions after reduction. + autoPtr<globalIndex> reduceRegions + ( + const label numLocalRegions, + const boolList& blockedFace, + labelList& cellRegion + ) const; + }; diff --git a/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.C b/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.C index d2c49eada057fb07139c386f56faedc4b2ffb9e8..c936e0b753ef3446fd3af109052c2fef15b206cc 100644 --- a/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.C +++ b/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.C @@ -239,7 +239,7 @@ Foam::labelList Foam::structuredRenumber::renumber deltaCalc.iterate(nLayers_); Info<< type() << " : did not visit " - << deltaCalc.getUnsetCells() + << deltaCalc.nUnvisitedCells() << " cells out of " << nTotalCells << "; using " << method_().type() << " renumbering for these" << endl;