diff --git a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H index b86e23c3c8..5a3386be77 100644 --- a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H +++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H @@ -25,10 +25,15 @@ labelList receptorNeigCell(mesh.nInternalFaces(), -1); PtrList meshParts(nZones); labelList nCellsPerZone(nZones, 0); - forAll(nCellsPerZone, zoneI) + // A mesh subset for each zone + forAll(meshParts, zonei) { - meshParts.set(zoneI, new fvMeshSubset(mesh)); - meshParts[zoneI].setLargeCellSubset(zoneID, zoneI); + meshParts.set + ( + zonei, + // Select cells where the zoneID == zonei + new fvMeshSubset(mesh, zonei, zoneID) + ); } for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++) diff --git a/applications/utilities/mesh/advanced/PDRMesh/PDRMesh.C b/applications/utilities/mesh/advanced/PDRMesh/PDRMesh.C index 549f6e58b0..e6547fb5b7 100644 --- a/applications/utilities/mesh/advanced/PDRMesh/PDRMesh.C +++ b/applications/utilities/mesh/advanced/PDRMesh/PDRMesh.C @@ -51,6 +51,7 @@ Description #include "fvMeshSubset.H" #include "argList.H" #include "cellSet.H" +#include "BitOps.H" #include "IOobjectList.H" #include "volFields.H" #include "mapPolyMesh.H" @@ -774,19 +775,19 @@ int main(int argc, char *argv[]) // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // - // Create mesh subsetting engine - fvMeshSubset subsetter(mesh); - - { - - cellSet blockedCells(mesh, blockedSetName); - - // invert - blockedCells.invert(mesh.nCells()); - - // Create subsetted mesh. - subsetter.setLargeCellSubset(blockedCells, defaultPatchi, true); - } + // Mesh subsetting engine + fvMeshSubset subsetter + ( + mesh, + BitSetOps::create + ( + mesh.nCells(), + cellSet(mesh, blockedSetName), // Blocked cells as labelHashSet + false // on=false: invert logic => retain the unblocked cells + ), + defaultPatchi, + true + ); // Subset wantedPatch. Note that might also include boundary faces diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C index 7301f82766..46668812b6 100644 --- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C +++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C @@ -570,17 +570,16 @@ labelList regionRenumber label celli = 0; - forAll(regionToCells, regionI) + forAll(regionToCells, regioni) { - Info<< " region " << regionI << " starts at " << celli << endl; + Info<< " region " << regioni << " starts at " << celli << endl; // Make sure no parallel comms - bool oldParRun = UPstream::parRun(); + const bool oldParRun = UPstream::parRun(); UPstream::parRun() = false; // Per region do a reordering. - fvMeshSubset subsetter(mesh); - subsetter.setLargeCellSubset(cellToRegion, regionI); + fvMeshSubset subsetter(mesh, regioni, cellToRegion); const fvMesh& subMesh = subsetter.subMesh(); diff --git a/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C b/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C index 565f048c5b..07fa8786ea 100644 --- a/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C +++ b/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation - \\/ M anipulation | Copyright (C) 2016-2017 OpenCFD Ltd. + \\/ M anipulation | Copyright (C) 2016-2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -45,6 +45,7 @@ Description #include "volFields.H" #include "topoDistanceData.H" #include "FaceCellWave.H" +#include "BitOps.H" #include "cellSet.H" #include "faceSet.H" #include "pointSet.H" @@ -311,7 +312,7 @@ void subsetTopoSets label nSet = 0; forAll(map, i) { - if (isSet[map[i]]) + if (isSet.test(map[i])) { ++nSet; } @@ -325,7 +326,7 @@ void subsetTopoSets TopoSet& subSet = subSets[i]; forAll(map, i) { - if (isSet[map[i]]) + if (isSet.test(map[i])) { subSet.insert(i); } @@ -371,7 +372,7 @@ int main(int argc, char *argv[]) #include "createNamedMesh.H" - const word setName = args[1]; + const word selectionName = args[1]; word meshInstance = mesh.pointsInstance(); word fieldsInstance = runTime.timeName(); @@ -389,30 +390,13 @@ int main(int argc, char *argv[]) } - Info<< "Reading cell set from " << setName << nl << endl; + Info<< "Reading cell set from " << selectionName << nl << endl; - // Create mesh subsetting engine - fvMeshSubset subsetter(mesh); - labelList exposedPatchIDs; + // Default exposed patch id + labelList exposedPatchIDs(one(), -1); - if (args.found("patch")) - { - const word patchName = args["patch"]; - - exposedPatchIDs = { mesh.boundaryMesh().findPatchID(patchName) }; - - if (exposedPatchIDs[0] == -1) - { - FatalErrorInFunction - << nl << "Valid patches are " << mesh.boundaryMesh().names() - << exit(FatalError); - } - - Info<< "Adding exposed internal faces to patch " << patchName - << nl << endl; - } - else if (args.found("patches")) + if (args.found("patches")) { const wordRes patchNames(args.readList("patches")); @@ -420,42 +404,77 @@ int main(int argc, char *argv[]) Info<< "Adding exposed internal faces to nearest of patches " << patchNames << nl << endl; + + if (exposedPatchIDs.empty()) + { + FatalErrorInFunction + << nl << "No patches matched. Patches: " + << mesh.boundaryMesh().names() + << exit(FatalError); + } + + } + else if (args.found("patch")) + { + const word patchName = args["patch"]; + + exposedPatchIDs.first() = mesh.boundaryMesh().findPatchID(patchName); + + Info<< "Adding exposed internal faces to patch " << patchName + << nl << endl; + + if (exposedPatchIDs.first() == -1) + { + FatalErrorInFunction + << nl << "No such patch. Patches: " + << mesh.boundaryMesh().names() + << exit(FatalError); + } } else { Info<< "Adding exposed internal faces to a patch called" << " \"oldInternalFaces\" (created if necessary)" << endl << endl; - exposedPatchIDs = { -1 }; } - cellSet currentSet(mesh, setName); + // Mesh subsetting engine + fvMeshSubset subsetter(mesh); - if (exposedPatchIDs.size() == 1) - { - subsetter.setLargeCellSubset(currentSet, exposedPatchIDs[0], true); - } - else - { - // Find per face the nearest patch - labelList nearestExposedPatch(nearestPatch(mesh, exposedPatchIDs)); + cellSet currentSet(mesh, selectionName); - labelList region(mesh.nCells(), 0); - forAllConstIter(cellSet, currentSet, iter) + { + bitSet selectedCells = BitSetOps::create(mesh.nCells(), currentSet); + + if (exposedPatchIDs.size() == 1) { - region[iter.key()] = 1; + // Single patch for exposed faces + subsetter.setCellSubset + ( + selectedCells, + exposedPatchIDs.first(), + true + ); } + else + { + // The nearest patch per face + labelList nearestExposedPatch(nearestPatch(mesh, exposedPatchIDs)); - labelList exposedFaces(subsetter.getExposedFaces(region, 1, true)); - subsetter.setLargeCellSubset - ( - region, - 1, - exposedFaces, - labelUIndList(nearestExposedPatch, exposedFaces)(), - true - ); + labelList exposedFaces + ( + subsetter.getExposedFaces(selectedCells, true) + ); + + subsetter.setCellSubset + ( + selectedCells, + exposedFaces, + labelUIndList(nearestExposedPatch, exposedFaces)(), + true + ); + } } diff --git a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C index f1b74976d2..07c493473b 100644 --- a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C +++ b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C @@ -872,30 +872,23 @@ autoPtr redistributeAndWrite // Find last non-processor patch. const polyBoundaryMesh& patches = mesh.boundaryMesh(); - label nonProcI = -1; + const label nonProcI = (patches.nNonProcessor() - 1); - forAll(patches, patchI) - { - if (isA(patches[patchI])) - { - break; - } - nonProcI++; - } - - if (nonProcI == -1) + if (nonProcI < 0) { FatalErrorInFunction << "Cannot find non-processor patch on processor " - << Pstream::myProcNo() << endl + << Pstream::myProcNo() << nl << " Current patches:" << patches.names() << abort(FatalError); } - // Subset 0 cells, no parallel comms. This is used to create - // zero-sized fields. - subsetterPtr.reset(new fvMeshSubset(mesh)); - subsetterPtr().setLargeCellSubset(labelHashSet(0), nonProcI, false); + // Subset 0 cells, no parallel comms. + // This is used to create zero-sized fields. + subsetterPtr.reset + ( + new fvMeshSubset(mesh, bitSet(), nonProcI, false) + ); } diff --git a/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoamMesh.C b/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoamMesh.C index 5970233049..7a29072240 100644 --- a/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoamMesh.C +++ b/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoamMesh.C @@ -423,8 +423,7 @@ void Foam::vtkPVFoam::convertMeshCellZones() if (!vtkgeom) { - fvMeshSubset subsetter(mesh); - subsetter.setLargeCellSubset(zMesh[zoneId]); + fvMeshSubset subsetter(mesh, zMesh[zoneId]); vtkgeom = vtuData.subset(subsetter, this->decomposePoly_); } @@ -490,8 +489,7 @@ void Foam::vtkPVFoam::convertMeshCellSets() if (!vtkgeom) { - fvMeshSubset subsetter(mesh); - subsetter.setLargeCellSubset(cellSet(mesh, partName)); + fvMeshSubset subsetter(mesh, cellSet(mesh, partName)); vtkgeom = vtuData.subset(subsetter, this->decomposePoly_); } diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C index e82e35a4fe..9d8ad5de9c 100644 --- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C +++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C @@ -2022,14 +2022,12 @@ Foam::autoPtr Foam::fvMeshDistribute::distribute //OPstream str(Pstream::commsTypes::blocking, recvProc); UOPstream str(recvProc, pBufs); - // Mesh subsetting engine - fvMeshSubset subsetter(mesh_); - - // Subset the cells of the current domain. - subsetter.setLargeCellSubset + // Mesh subsetting engine - subset the cells of the current domain. + fvMeshSubset subsetter ( - distribution, + mesh_, recvProc, + distribution, oldInternalPatchi, // oldInternalFaces patch false // no parallel sync ); diff --git a/src/dynamicMesh/fvMeshSubset/fvMeshSubset.C b/src/dynamicMesh/fvMeshSubset/fvMeshSubset.C index 48b9278069..239f9e1704 100644 --- a/src/dynamicMesh/fvMeshSubset/fvMeshSubset.C +++ b/src/dynamicMesh/fvMeshSubset/fvMeshSubset.C @@ -21,15 +21,12 @@ License You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see . -Description - Post-processing mesh subset tool. Given the original mesh and the - list of selected cells, it creates the mesh consisting only of the - desired cells, with the mapping list for points, faces, and cells. - \*---------------------------------------------------------------------------*/ #include "fvMeshSubset.H" #include "boolList.H" +#include "BitOps.H" +#include "pointIndList.H" #include "Pstream.H" #include "emptyPolyPatch.H" #include "demandDrivenData.H" @@ -38,6 +35,27 @@ Description #include "polyTopoChange.H" #include "mapPolyMesh.H" +// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // + +namespace +{ + +// Mark faces/points (with 0) in labelList +inline void markUsed +( + const Foam::labelUList& locations, + Foam::labelList& map +) +{ + for (auto idx : locations) + { + map[idx] = 0; + } +} + +} // End anonymous namespace + + // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // bool Foam::fvMeshSubset::checkCellSubset() const @@ -56,15 +74,37 @@ bool Foam::fvMeshSubset::checkCellSubset() const } -void Foam::fvMeshSubset::markPoints -( - const labelUList& curPoints, - labelList& pointMap -) +void Foam::fvMeshSubset::calcFaceFlipMap() const { - for (const label pointi : curPoints) + const labelList& subToBaseFace = faceMap(); + const labelList& subToBaseCell = cellMap(); + + faceFlipMapPtr_.clear(); + faceFlipMapPtr_.reset(new labelList(subToBaseFace.size())); + labelList& faceFlipMap = *faceFlipMapPtr_; + + // Only exposed internal faces might be flipped (since we don't do + // any cell renumbering, just compacting) + const label subInt = subMesh().nInternalFaces(); + + const labelList& subOwn = subMesh().faceOwner(); + const labelList& own = baseMesh_.faceOwner(); + + for (label subFaceI = 0; subFaceI < subInt; ++subFaceI) { - pointMap[pointi] = 0; + faceFlipMap[subFaceI] = subToBaseFace[subFaceI]+1; + } + for (label subFaceI = subInt; subFaceI < subOwn.size(); ++subFaceI) + { + const label faceI = subToBaseFace[subFaceI]; + if (subToBaseCell[subOwn[subFaceI]] == own[faceI]) + { + faceFlipMap[subFaceI] = faceI+1; + } + else + { + faceFlipMap[subFaceI] = -faceI-1; + } } } @@ -87,10 +127,8 @@ void Foam::fvMeshSubset::doCoupledPatches PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking); // Send face usage across processor patches - forAll(oldPatches, oldPatchi) + for (const polyPatch& pp : oldPatches) { - const polyPatch& pp = oldPatches[oldPatchi]; - if (isA(pp)) { const processorPolyPatch& procPatch = @@ -113,10 +151,8 @@ void Foam::fvMeshSubset::doCoupledPatches pBufs.finishedSends(); // Receive face usage count and check for faces that become uncoupled. - forAll(oldPatches, oldPatchi) + for (const polyPatch& pp : oldPatches) { - const polyPatch& pp = oldPatches[oldPatchi]; - if (isA(pp)) { const processorPolyPatch& procPatch = @@ -145,7 +181,7 @@ void Foam::fvMeshSubset::doCoupledPatches // Face's neighbour is no longer there. Mark face // off as coupled nCellsUsingFace[pp.start()+i] = 3; - nUncoupled++; + ++nUncoupled; } } } @@ -154,10 +190,8 @@ void Foam::fvMeshSubset::doCoupledPatches } // Do same for cyclics. - forAll(oldPatches, oldPatchi) + for (const polyPatch& pp : oldPatches) { - const polyPatch& pp = oldPatches[oldPatchi]; - if (isA(pp)) { const cyclicPolyPatch& cycPatch = @@ -177,7 +211,7 @@ void Foam::fvMeshSubset::doCoupledPatches ) { nCellsUsingFace[thisFacei] = 3; - nUncoupled++; + ++nUncoupled; } } } @@ -197,862 +231,14 @@ void Foam::fvMeshSubset::doCoupledPatches } -Foam::labelList Foam::fvMeshSubset::subset +void Foam::fvMeshSubset::removeCellsImpl ( - const label nElems, - const labelList& selectedElements, - const labelList& subsetMap -) -{ - // Mark selected elements. - boolList selected(nElems, false); - forAll(selectedElements, i) - { - selected[selectedElements[i]] = true; - } - - // Count subset of selected elements - label n = 0; - forAll(subsetMap, i) - { - if (selected[subsetMap[i]]) - { - n++; - } - } - - // Collect selected elements - labelList subsettedElements(n); - n = 0; - - forAll(subsetMap, i) - { - if (selected[subsetMap[i]]) - { - subsettedElements[n++] = i; - } - } - - return subsettedElements; -} - - -void Foam::fvMeshSubset::subsetZones() -{ - // Keep all zones, even if zero size. - - const pointZoneMesh& pointZones = baseMesh().pointZones(); - - // PointZones - List pZonePtrs(pointZones.size()); - - forAll(pointZones, i) - { - const pointZone& pz = pointZones[i]; - - pZonePtrs[i] = new pointZone - ( - pz.name(), - subset(baseMesh().nPoints(), pz, pointMap()), - i, - fvMeshSubsetPtr_().pointZones() - ); - } - - - // FaceZones - - const faceZoneMesh& faceZones = baseMesh().faceZones(); - - - // Do we need to remove zones where the side we're interested in - // no longer exists? Guess not. - List fZonePtrs(faceZones.size()); - - forAll(faceZones, i) - { - const faceZone& fz = faceZones[i]; - - // Expand faceZone to full mesh - // +1 : part of faceZone, flipped - // -1 : ,, , unflipped - // 0 : not part of faceZone - labelList zone(baseMesh().nFaces(), 0); - forAll(fz, j) - { - if (fz.flipMap()[j]) - { - zone[fz[j]] = 1; - } - else - { - zone[fz[j]] = -1; - } - } - - // Select faces - label nSub = 0; - forAll(faceMap(), j) - { - if (zone[faceMap()[j]] != 0) - { - nSub++; - } - } - labelList subAddressing(nSub); - boolList subFlipStatus(nSub); - nSub = 0; - forAll(faceMap(), subFacei) - { - label meshFacei = faceMap()[subFacei]; - if (zone[meshFacei] != 0) - { - subAddressing[nSub] = subFacei; - label subOwner = subMesh().faceOwner()[subFacei]; - label baseOwner = baseMesh().faceOwner()[meshFacei]; - // If subowner is the same cell as the base keep the flip status - bool sameOwner = (cellMap()[subOwner] == baseOwner); - bool flip = (zone[meshFacei] == 1); - subFlipStatus[nSub] = (sameOwner == flip); - - nSub++; - } - } - - fZonePtrs[i] = new faceZone - ( - fz.name(), - subAddressing, - subFlipStatus, - i, - fvMeshSubsetPtr_().faceZones() - ); - } - - - const cellZoneMesh& cellZones = baseMesh().cellZones(); - - List cZonePtrs(cellZones.size()); - - forAll(cellZones, i) - { - const cellZone& cz = cellZones[i]; - - cZonePtrs[i] = new cellZone - ( - cz.name(), - subset(baseMesh().nCells(), cz, cellMap()), - i, - fvMeshSubsetPtr_().cellZones() - ); - } - - - // Add the zones - fvMeshSubsetPtr_().addZones(pZonePtrs, fZonePtrs, cZonePtrs); -} - - -Foam::labelList Foam::fvMeshSubset::getCellsToRemove -( - const labelList& region, - const label currentRegion -) const -{ - // Count - label nKeep = 0; - forAll(region, cellI) - { - if (region[cellI] == currentRegion) - { - nKeep++; - } - } - - // Collect cells to remove - label nRemove = baseMesh().nCells() - nKeep; - labelList cellsToRemove(nRemove); - - nRemove = 0; - forAll(region, cellI) - { - if (region[cellI] != currentRegion) - { - cellsToRemove[nRemove++] = cellI; - } - } - - return cellsToRemove; -} - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::fvMeshSubset::fvMeshSubset(const fvMesh& baseMesh) -: - baseMesh_(baseMesh), - fvMeshSubsetPtr_(nullptr), - pointMap_(0), - faceMap_(0), - cellMap_(0), - patchMap_(0), - faceFlipMapPtr_() -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -void Foam::fvMeshSubset::setLargeCellSubset -( - const labelList& region, - const label currentRegion, - const label patchID, - const bool syncPar -) -{ - const cellList& oldCells = baseMesh().cells(); - const faceList& oldFaces = baseMesh().faces(); - const pointField& oldPoints = baseMesh().points(); - const labelList& oldOwner = baseMesh().faceOwner(); - const labelList& oldNeighbour = baseMesh().faceNeighbour(); - const polyBoundaryMesh& oldPatches = baseMesh().boundaryMesh(); - const label oldNInternalFaces = baseMesh().nInternalFaces(); - - // Initial checks - - if (region.size() != oldCells.size()) - { - FatalErrorInFunction - << "Size of region " << region.size() - << " is not equal to number of cells in mesh " << oldCells.size() - << abort(FatalError); - } - - // Initial check on patches before doing anything time consuming. - - label wantedPatchID = patchID; - - if (wantedPatchID == -1) - { - // No explicit patch specified. Put in oldInternalFaces patch. - // Check if patch with this name already exists. - wantedPatchID = oldPatches.findPatchID("oldInternalFaces"); - } - else if (wantedPatchID < 0 || wantedPatchID >= oldPatches.size()) - { - FatalErrorInFunction - << "Non-existing patch index " << wantedPatchID << endl - << "Should be between 0 and " << oldPatches.size()-1 - << abort(FatalError); - } - - // Clear demand driven data - faceFlipMapPtr_.clear(); - - // Get the cells for the current region - sorted in ascending order - cellMap_.setSize(oldCells.size()); - - { - label n = 0; - forAll(oldCells, oldCelli) - { - if (region[oldCelli] == currentRegion) - { - cellMap_[n++] = oldCelli; - } - } - cellMap_.setSize(n); - } - - - // Mark all used faces. Count number of cells using them - // 0: face not used anymore - // 1: face used by one cell, face becomes/stays boundary face - // 2: face still used and remains internal face - // 3: face coupled and used by one cell only (so should become normal, - // non-coupled patch face) - // - // Note that this is not really necessary - but means we can size things - // correctly. Also makes handling coupled faces much easier. - - labelList nCellsUsingFace(oldFaces.size(), 0); - - label nFacesInSet = 0; - forAll(oldFaces, oldFacei) - { - bool faceUsed = false; - - if (region[oldOwner[oldFacei]] == currentRegion) - { - nCellsUsingFace[oldFacei]++; - faceUsed = true; - } - - if - ( - baseMesh().isInternalFace(oldFacei) - && (region[oldNeighbour[oldFacei]] == currentRegion) - ) - { - nCellsUsingFace[oldFacei]++; - faceUsed = true; - } - - if (faceUsed) - { - nFacesInSet++; - } - } - faceMap_.setSize(nFacesInSet); - - // Handle coupled faces. Modifies patch faces to be uncoupled to 3. - doCoupledPatches(syncPar, nCellsUsingFace); - - - // See which patch to use for exposed internal faces. - label oldInternalPatchID = 0; - - // Insert faces before which patch - label nextPatchID = oldPatches.size(); - - // old to new patches - labelList globalPatchMap(oldPatches.size()); - - // New patch size - label nbSize = oldPatches.size(); - - if (wantedPatchID == -1) - { - // Create 'oldInternalFaces' patch at the end (or before - // processorPatches) - // and put all exposed internal faces in there. - - forAll(oldPatches, patchi) - { - if (isA(oldPatches[patchi])) - { - nextPatchID = patchi; - break; - } - oldInternalPatchID++; - } - - nbSize++; - - // adapt old to new patches for inserted patch - for (label oldPatchi = 0; oldPatchi < nextPatchID; oldPatchi++) - { - globalPatchMap[oldPatchi] = oldPatchi; - } - for - ( - label oldPatchi = nextPatchID; - oldPatchi < oldPatches.size(); - oldPatchi++ - ) - { - globalPatchMap[oldPatchi] = oldPatchi+1; - } - } - else - { - oldInternalPatchID = wantedPatchID; - nextPatchID = wantedPatchID+1; - - // old to new patches - globalPatchMap = identity(oldPatches.size()); - } - - labelList boundaryPatchSizes(nbSize, 0); - - - // Make a global-to-local point map - labelList globalPointMap(oldPoints.size(), -1); - - labelList globalFaceMap(oldFaces.size(), -1); - label facei = 0; - - // 1. Pick up all preserved internal faces. - for (label oldFacei = 0; oldFacei < oldNInternalFaces; oldFacei++) - { - if (nCellsUsingFace[oldFacei] == 2) - { - globalFaceMap[oldFacei] = facei; - faceMap_[facei++] = oldFacei; - - // Mark all points from the face - markPoints(oldFaces[oldFacei], globalPointMap); - } - } - - // These are all the internal faces in the mesh. - label nInternalFaces = facei; - - // 2. Boundary faces up to where we want to insert old internal faces - for - ( - label oldPatchi = 0; - oldPatchi < oldPatches.size() - && oldPatchi < nextPatchID; - oldPatchi++ - ) - { - const polyPatch& oldPatch = oldPatches[oldPatchi]; - - label oldFacei = oldPatch.start(); - - forAll(oldPatch, i) - { - if (nCellsUsingFace[oldFacei] == 1) - { - // Boundary face is kept. - - // Mark face and increment number of points in set - globalFaceMap[oldFacei] = facei; - faceMap_[facei++] = oldFacei; - - // Mark all points from the face - markPoints(oldFaces[oldFacei], globalPointMap); - - // Increment number of patch faces - boundaryPatchSizes[globalPatchMap[oldPatchi]]++; - } - oldFacei++; - } - } - - // 3a. old internal faces that have become exposed. - for (label oldFacei = 0; oldFacei < oldNInternalFaces; oldFacei++) - { - if (nCellsUsingFace[oldFacei] == 1) - { - globalFaceMap[oldFacei] = facei; - faceMap_[facei++] = oldFacei; - - // Mark all points from the face - markPoints(oldFaces[oldFacei], globalPointMap); - - // Increment number of patch faces - boundaryPatchSizes[oldInternalPatchID]++; - } - } - - // 3b. coupled patch faces that have become uncoupled. - for - ( - label oldFacei = oldNInternalFaces; - oldFacei < oldFaces.size(); - oldFacei++ - ) - { - if (nCellsUsingFace[oldFacei] == 3) - { - globalFaceMap[oldFacei] = facei; - faceMap_[facei++] = oldFacei; - - // Mark all points from the face - markPoints(oldFaces[oldFacei], globalPointMap); - - // Increment number of patch faces - boundaryPatchSizes[oldInternalPatchID]++; - } - } - - // 4. Remaining boundary faces - for - ( - label oldPatchi = nextPatchID; - oldPatchi < oldPatches.size(); - oldPatchi++ - ) - { - const polyPatch& oldPatch = oldPatches[oldPatchi]; - - label oldFacei = oldPatch.start(); - - forAll(oldPatch, i) - { - if (nCellsUsingFace[oldFacei] == 1) - { - // Boundary face is kept. - - // Mark face and increment number of points in set - globalFaceMap[oldFacei] = facei; - faceMap_[facei++] = oldFacei; - - // Mark all points from the face - markPoints(oldFaces[oldFacei], globalPointMap); - - // Increment number of patch faces - boundaryPatchSizes[globalPatchMap[oldPatchi]]++; - } - oldFacei++; - } - } - - if (facei != nFacesInSet) - { - FatalErrorInFunction - << "Problem" << abort(FatalError); - } - - - // Grab the points map - label nPointsInSet = 0; - - forAll(globalPointMap, pointi) - { - if (globalPointMap[pointi] != -1) - { - nPointsInSet++; - } - } - pointMap_.setSize(nPointsInSet); - - nPointsInSet = 0; - - forAll(globalPointMap, pointi) - { - if (globalPointMap[pointi] != -1) - { - pointMap_[nPointsInSet] = pointi; - globalPointMap[pointi] = nPointsInSet; - nPointsInSet++; - } - } - - //Pout<< "Number of cells in new mesh : " << cellMap_.size() << endl; - //Pout<< "Number of faces in new mesh : " << faceMap_.size() << endl; - //Pout<< "Number of points in new mesh: " << pointMap_.size() << endl; - - // Make a new mesh - pointField newPoints(pointMap_.size()); - - label nNewPoints = 0; - - forAll(pointMap_, pointi) - { - newPoints[nNewPoints] = oldPoints[pointMap_[pointi]]; - nNewPoints++; - } - - faceList newFaces(faceMap_.size()); - - label nNewFaces = 0; - - // Make internal faces - for (label facei = 0; facei < nInternalFaces; facei++) - { - const face& oldF = oldFaces[faceMap_[facei]]; - - face newF(oldF.size()); - - forAll(newF, i) - { - newF[i] = globalPointMap[oldF[i]]; - } - - newFaces[nNewFaces] = newF; - nNewFaces++; - } - - - // Make boundary faces. (different from internal since might need to be - // flipped) - for (label facei = nInternalFaces; facei < faceMap_.size(); facei++) - { - label oldFacei = faceMap_[facei]; - - face oldF = oldFaces[oldFacei]; - - // Turn the faces as necessary to point outwards - if (baseMesh().isInternalFace(oldFacei)) - { - // Was internal face. Possibly turned the wrong way round - if - ( - region[oldOwner[oldFacei]] != currentRegion - && region[oldNeighbour[oldFacei]] == currentRegion - ) - { - oldF = oldFaces[oldFacei].reverseFace(); - } - } - - // Relabel vertices of the (possibly turned) face. - face newF(oldF.size()); - - forAll(newF, i) - { - newF[i] = globalPointMap[oldF[i]]; - } - - newFaces[nNewFaces] = newF; - nNewFaces++; - } - - - // Create cells - cellList newCells(cellMap_.size()); - - label nNewCells = 0; - - forAll(cellMap_, celli) - { - const labelList& oldC = oldCells[cellMap_[celli]]; - - labelList newC(oldC.size()); - - forAll(newC, i) - { - newC[i] = globalFaceMap[oldC[i]]; - } - - newCells[nNewCells] = cell(newC); - nNewCells++; - } - - - // Delete any old mesh first - fvMeshSubsetPtr_.clear(); - - // Make a new mesh - // Note that mesh gets registered with same name as original mesh. This is - // not proper but cannot be avoided since otherwise surfaceInterpolation - // cannot find its fvSchemes (it will try to read e.g. - // system/region0SubSet/fvSchemes) - fvMeshSubsetPtr_ = autoPtr::New - ( - IOobject - ( - baseMesh().name(), - baseMesh().time().timeName(), - baseMesh().time(), - IOobject::NO_READ, - IOobject::NO_WRITE - ), - std::move(newPoints), - std::move(newFaces), - std::move(newCells), - syncPar // parallel synchronisation - ); - - // Add old patches - List newBoundary(nbSize); - patchMap_.setSize(nbSize); - label nNewPatches = 0; - label patchStart = nInternalFaces; - - - // For parallel: only remove patch if none of the processors has it. - // This only gets done for patches before the one being inserted - // (so patches < nextPatchID) - - // Get sum of patch sizes. Zero if patch can be deleted. - labelList globalPatchSizes(boundaryPatchSizes); - globalPatchSizes.setSize(nextPatchID); - - if (syncPar && Pstream::parRun()) - { - // Get patch names (up to nextPatchID) - List patchNames(Pstream::nProcs()); - patchNames[Pstream::myProcNo()] = oldPatches.names(); - patchNames[Pstream::myProcNo()].setSize(nextPatchID); - Pstream::gatherList(patchNames); - Pstream::scatterList(patchNames); - - // Get patch sizes (up to nextPatchID). - // Note that up to nextPatchID the globalPatchMap is an identity so - // no need to index through that. - Pstream::listCombineGather(globalPatchSizes, plusEqOp