48 bool Foam::fvMeshSubset::checkCellSubset()
const 50 if (fvMeshSubsetPtr_.empty())
53 <<
"void setCellSubset(const labelHashSet& cellsToSubset)" <<
endl 54 <<
"before attempting to access subset data" 66 void Foam::fvMeshSubset::markPoints
75 pointMap.insert(curPoints[pointi], 0);
80 void Foam::fvMeshSubset::markPoints
88 pointMap[curPoints[pointi]] = 0;
93 void Foam::fvMeshSubset::doCoupledPatches
102 const polyBoundaryMesh& oldPatches = baseMesh().boundaryMesh();
104 label nUncoupled = 0;
111 forAll(oldPatches, oldPatchi)
113 const polyPatch& pp = oldPatches[oldPatchi];
115 if (isA<processorPolyPatch>(pp))
117 const processorPolyPatch& procPatch =
118 refCast<const processorPolyPatch>(pp);
120 UOPstream toNeighbour(procPatch.neighbProcNo(), pBufs);
123 << SubList<label>(nCellsUsingFace, pp.size(), pp.start());
127 pBufs.finishedSends();
130 forAll(oldPatches, oldPatchi)
132 const polyPatch& pp = oldPatches[oldPatchi];
134 if (isA<processorPolyPatch>(pp))
136 const processorPolyPatch& procPatch =
137 refCast<const processorPolyPatch>(pp);
139 UIPstream fromNeighbour(procPatch.neighbProcNo(), pBufs);
141 labelList nbrCellsUsingFace(fromNeighbour);
149 nCellsUsingFace[pp.start()+i] == 1
150 && nbrCellsUsingFace[i] == 0
155 nCellsUsingFace[pp.start()+i] = 3;
164 forAll(oldPatches, oldPatchi)
166 const polyPatch& pp = oldPatches[oldPatchi];
168 if (isA<cyclicPolyPatch>(pp))
170 const cyclicPolyPatch& cycPatch =
171 refCast<const cyclicPolyPatch>(pp);
175 label thisFacei = cycPatch.start() + i;
176 label otherFacei = cycPatch.transformGlobalFace(thisFacei);
180 nCellsUsingFace[thisFacei] == 1
181 && nCellsUsingFace[otherFacei] == 0
184 nCellsUsingFace[thisFacei] = 3;
193 reduce(nUncoupled, sumOp<label>());
198 Info<<
"Uncoupled " << nUncoupled <<
" faces on coupled patches. " 199 <<
"(processorPolyPatch, cyclicPolyPatch)" <<
endl;
213 forAll(selectedElements, i)
215 selected[selectedElements[i]] =
true;
222 if (selected[subsetMap[i]])
234 if (selected[subsetMap[i]])
236 subsettedElements[n++] = i;
240 return subsettedElements;
244 void Foam::fvMeshSubset::subsetZones()
251 List<pointZone*> pZonePtrs(pointZones.size());
255 const pointZone& pz = pointZones[i];
257 pZonePtrs[i] =
new pointZone
262 fvMeshSubsetPtr_().pointZones()
274 List<faceZone*> fZonePtrs(faceZones.size());
278 const faceZone& fz = faceZones[i];
312 if (zone[meshFacei] != 0)
314 subAddressing[nSub] = subFacei;
315 label subOwner = subMesh().faceOwner()[subFacei];
316 label baseOwner = baseMesh().faceOwner()[meshFacei];
318 bool sameOwner = (cellMap()[subOwner] == baseOwner);
319 bool flip = (zone[meshFacei] == 1);
320 subFlipStatus[nSub] = (sameOwner == flip);
326 fZonePtrs[i] =
new faceZone
332 fvMeshSubsetPtr_().faceZones()
339 List<cellZone*> cZonePtrs(cellZones.size());
343 const cellZone& cz = cellZones[i];
345 cZonePtrs[i] =
new cellZone
348 subset(baseMesh().nCells(), cz, cellMap()),
350 fvMeshSubsetPtr_().cellZones()
356 fvMeshSubsetPtr_().addZones(pZonePtrs, fZonePtrs, cZonePtrs);
363 const label currentRegion
370 if (region[cellI] == currentRegion)
377 label nRemove = baseMesh().nCells() - nKeep;
383 if (region[cellI] != currentRegion)
385 cellsToRemove[nRemove++] = cellI;
389 return cellsToRemove;
395 Foam::fvMeshSubset::fvMeshSubset(
const fvMesh& baseMesh)
398 fvMeshSubsetPtr_(nullptr),
423 label wantedPatchID = patchID;
425 if (wantedPatchID == -1)
429 wantedPatchID = oldPatches.
findPatchID(
"oldInternalFaces");
431 else if (wantedPatchID < 0 || wantedPatchID >= oldPatches.
size())
434 <<
"Non-existing patch index " << wantedPatchID <<
endl 435 <<
"Should be between 0 and " << oldPatches.
size()-1
441 faceFlipMapPtr_.clear();
444 cellMap_ = globalCellMap.
toc();
450 const label avgNFacesPerCell = 6;
451 const label avgNPointsPerFace = 4;
458 Map<label> facesToSubset(avgNFacesPerCell*nCellsInSet);
463 const labelList& curFaces = oldCells[cellMap_[celli]];
467 if (!facesToSubset.
found(curFaces[facei]))
469 facesToSubset.
insert(curFaces[facei], 1);
473 facesToSubset[curFaces[facei]]++;
482 Map<label> globalPointMap(avgNPointsPerFace*facesToSubset.
size());
494 if (facesToSubset[facesToc[facei]] == 2)
497 faceMap_[globalFaceMap.size()] = facesToc[facei];
498 globalFaceMap.
insert(facesToc[facei], globalFaceMap.
size());
501 markPoints(oldFaces[facesToc[facei]], globalPointMap);
506 label nInternalFaces = globalFaceMap.size();
511 if (wantedPatchID != -1)
513 oldPatchStart = oldPatches[wantedPatchID].start();
520 for (; facei< facesToc.
size(); facei++)
522 if (facesToc[facei] >= oldPatchStart)
529 && facesToSubset[facesToc[facei]] == 1
533 faceMap_[globalFaceMap.size()] = facesToc[facei];
534 globalFaceMap.insert(facesToc[facei], globalFaceMap.size());
537 markPoints(oldFaces[facesToc[facei]], globalPointMap);
542 forAll(facesToc, intFacei)
546 baseMesh().isInternalFace(facesToc[intFacei])
547 && facesToSubset[facesToc[intFacei]] == 1
551 faceMap_[globalFaceMap.size()] = facesToc[intFacei];
552 globalFaceMap.
insert(facesToc[intFacei], globalFaceMap.
size());
555 markPoints(oldFaces[facesToc[intFacei]], globalPointMap);
560 for (; facei< facesToc.
size(); facei++)
565 && facesToSubset[facesToc[facei]] == 1
569 faceMap_[globalFaceMap.size()] = facesToc[facei];
570 globalFaceMap.insert(facesToc[facei], globalFaceMap.size());
573 markPoints(oldFaces[facesToc[facei]], globalPointMap);
580 pointMap_ = globalPointMap.toc();
585 globalPointMap[pointMap_[pointi]] = pointi;
588 Pout<<
"Number of cells in new mesh: " << nCellsInSet <<
endl;
589 Pout<<
"Number of faces in new mesh: " << globalFaceMap.size() <<
endl;
590 Pout<<
"Number of points in new mesh: " << globalPointMap.size() <<
endl;
595 label nNewPoints = 0;
599 newPoints[nNewPoints] = oldPoints[pointMap_[pointi]];
603 faceList newFaces(globalFaceMap.size());
608 for (
label facei = 0; facei < nInternalFaces; facei++)
610 const face& oldF = oldFaces[faceMap_[facei]];
616 newF[i] = globalPointMap[oldF[i]];
619 newFaces[nNewFaces] = newF;
626 label oldInternalPatchID = -1;
628 if (wantedPatchID == -1)
632 oldInternalPatchID = nbSize;
638 oldInternalPatchID = wantedPatchID;
648 for (
label facei = nInternalFaces; facei < faceMap_.
size(); facei++)
650 label oldFacei = faceMap_[facei];
652 face oldF = oldFaces[oldFacei];
660 !globalCellMap.
found(oldOwner[oldFacei])
661 && globalCellMap.
found(oldNeighbour[oldFacei])
664 oldF = oldFaces[oldFacei].reverseFace();
668 boundaryPatchSizes[oldInternalPatchID]++;
676 boundaryPatchSizes[patchOfFace]++;
683 newF[i] = globalPointMap[oldF[i]];
686 newFaces[nNewFaces] = newF;
699 const labelList& oldC = oldCells[cellMap_[celli]];
705 newC[i] = globalFaceMap[oldC[i]];
708 newCells[nNewCells] =
cell(newC);
714 fvMeshSubsetPtr_.clear();
716 fvMeshSubsetPtr_.reset
738 label nNewPatches = 0;
739 label patchStart = nInternalFaces;
743 if (boundaryPatchSizes[
patchi] > 0)
746 newBoundary[nNewPatches] = oldPatches[
patchi].
clone 750 boundaryPatchSizes[
patchi],
754 patchStart += boundaryPatchSizes[
patchi];
755 patchMap_[nNewPatches] =
patchi;
760 if (wantedPatchID == -1)
763 if (boundaryPatchSizes[oldInternalPatchID] > 0)
768 boundaryPatchSizes[oldInternalPatchID],
772 emptyPolyPatch::typeName
777 patchMap_[nNewPatches] = -1;
783 newBoundary.
setSize(nNewPatches);
784 patchMap_.
setSize(nNewPatches);
787 fvMeshSubsetPtr_().addFvPatches(newBoundary);
797 const label currentRegion,
812 if (region.
size() != oldCells.
size())
815 <<
"Size of region " << region.
size()
816 <<
" is not equal to number of cells in mesh " << oldCells.
size()
821 label wantedPatchID = patchID;
823 if (wantedPatchID == -1)
827 wantedPatchID = oldPatches.
findPatchID(
"oldInternalFaces");
829 else if (wantedPatchID < 0 || wantedPatchID >= oldPatches.
size())
832 <<
"Non-existing patch index " << wantedPatchID <<
endl 833 <<
"Should be between 0 and " << oldPatches.
size()-1
838 faceFlipMapPtr_.clear();
842 label nCellsInSet = 0;
846 if (region[oldCelli] == currentRegion)
848 cellMap_[nCellsInSet++] = oldCelli;
866 label nFacesInSet = 0;
867 forAll(oldFaces, oldFacei)
869 bool faceUsed =
false;
871 if (region[oldOwner[oldFacei]] == currentRegion)
873 nCellsUsingFace[oldFacei]++;
880 && (region[oldNeighbour[oldFacei]] == currentRegion)
883 nCellsUsingFace[oldFacei]++;
895 doCoupledPatches(syncPar, nCellsUsingFace);
899 label oldInternalPatchID = 0;
910 if (wantedPatchID == -1)
918 if (isA<processorPolyPatch>(oldPatches[
patchi]))
923 oldInternalPatchID++;
929 for (
label oldPatchi = 0; oldPatchi < nextPatchID; oldPatchi++)
931 globalPatchMap[oldPatchi] = oldPatchi;
935 label oldPatchi = nextPatchID;
936 oldPatchi < oldPatches.
size();
940 globalPatchMap[oldPatchi] = oldPatchi+1;
945 oldInternalPatchID = wantedPatchID;
946 nextPatchID = wantedPatchID+1;
962 for (
label oldFacei = 0; oldFacei < oldNInternalFaces; oldFacei++)
964 if (nCellsUsingFace[oldFacei] == 2)
966 globalFaceMap[oldFacei] = facei;
967 faceMap_[facei++] = oldFacei;
970 markPoints(oldFaces[oldFacei], globalPointMap);
975 label nInternalFaces = facei;
981 oldPatchi < oldPatches.
size()
982 && oldPatchi < nextPatchID;
986 const polyPatch& oldPatch = oldPatches[oldPatchi];
992 if (nCellsUsingFace[oldFacei] == 1)
997 globalFaceMap[oldFacei] = facei;
998 faceMap_[facei++] = oldFacei;
1001 markPoints(oldFaces[oldFacei], globalPointMap);
1004 boundaryPatchSizes[globalPatchMap[oldPatchi]]++;
1011 for (
label oldFacei = 0; oldFacei < oldNInternalFaces; oldFacei++)
1013 if (nCellsUsingFace[oldFacei] == 1)
1015 globalFaceMap[oldFacei] = facei;
1016 faceMap_[facei++] = oldFacei;
1019 markPoints(oldFaces[oldFacei], globalPointMap);
1022 boundaryPatchSizes[oldInternalPatchID]++;
1029 label oldFacei = oldNInternalFaces;
1030 oldFacei < oldFaces.
size();
1034 if (nCellsUsingFace[oldFacei] == 3)
1036 globalFaceMap[oldFacei] = facei;
1037 faceMap_[facei++] = oldFacei;
1040 markPoints(oldFaces[oldFacei], globalPointMap);
1043 boundaryPatchSizes[oldInternalPatchID]++;
1050 label oldPatchi = nextPatchID;
1051 oldPatchi < oldPatches.
size();
1055 const polyPatch& oldPatch = oldPatches[oldPatchi];
1061 if (nCellsUsingFace[oldFacei] == 1)
1066 globalFaceMap[oldFacei] = facei;
1067 faceMap_[facei++] = oldFacei;
1070 markPoints(oldFaces[oldFacei], globalPointMap);
1073 boundaryPatchSizes[globalPatchMap[oldPatchi]]++;
1079 if (facei != nFacesInSet)
1087 label nPointsInSet = 0;
1089 forAll(globalPointMap, pointi)
1091 if (globalPointMap[pointi] != -1)
1096 pointMap_.
setSize(nPointsInSet);
1100 forAll(globalPointMap, pointi)
1102 if (globalPointMap[pointi] != -1)
1104 pointMap_[nPointsInSet] = pointi;
1105 globalPointMap[pointi] = nPointsInSet;
1117 label nNewPoints = 0;
1119 forAll(pointMap_, pointi)
1121 newPoints[nNewPoints] = oldPoints[pointMap_[pointi]];
1127 label nNewFaces = 0;
1130 for (
label facei = 0; facei < nInternalFaces; facei++)
1132 const face& oldF = oldFaces[faceMap_[facei]];
1138 newF[i] = globalPointMap[oldF[i]];
1141 newFaces[nNewFaces] = newF;
1148 for (
label facei = nInternalFaces; facei < faceMap_.
size(); facei++)
1150 label oldFacei = faceMap_[facei];
1152 face oldF = oldFaces[oldFacei];
1160 region[oldOwner[oldFacei]] != currentRegion
1161 && region[oldNeighbour[oldFacei]] == currentRegion
1164 oldF = oldFaces[oldFacei].reverseFace();
1173 newF[i] = globalPointMap[oldF[i]];
1176 newFaces[nNewFaces] = newF;
1185 label nNewCells = 0;
1189 const labelList& oldC = oldCells[cellMap_[celli]];
1195 newC[i] = globalFaceMap[oldC[i]];
1198 newCells[nNewCells] =
cell(newC);
1204 fvMeshSubsetPtr_.clear();
1212 fvMeshSubsetPtr_.reset
1234 label nNewPatches = 0;
1235 label patchStart = nInternalFaces;
1243 labelList globalPatchSizes(boundaryPatchSizes);
1244 globalPatchSizes.
setSize(nextPatchID);
1264 bool samePatches =
true;
1266 for (
label proci = 1; proci < patchNames.
size(); proci++)
1268 if (patchNames[proci] != patchNames[0])
1270 samePatches =
false;
1288 label oldPatchi = 0;
1289 oldPatchi < oldPatches.
size()
1290 && oldPatchi < nextPatchID;
1294 label newSize = boundaryPatchSizes[globalPatchMap[oldPatchi]];
1297 newBoundary[nNewPatches] = oldPatches[oldPatchi].
clone 1305 patchStart += newSize;
1306 patchMap_[nNewPatches] = oldPatchi;
1312 if (wantedPatchID == -1)
1314 label oldInternalSize = boundaryPatchSizes[oldInternalPatchID];
1322 if (oldInternalSize > 0)
1327 boundaryPatchSizes[oldInternalPatchID],
1331 emptyPolyPatch::typeName
1339 patchStart += boundaryPatchSizes[oldInternalPatchID];
1340 patchMap_[nNewPatches] = -1;
1349 label oldPatchi = nextPatchID;
1350 oldPatchi < oldPatches.
size();
1354 label newSize = boundaryPatchSizes[globalPatchMap[oldPatchi]];
1357 newBoundary[nNewPatches] = oldPatches[oldPatchi].
clone 1368 patchStart += newSize;
1369 patchMap_[nNewPatches] = oldPatchi;
1375 newBoundary.
setSize(nNewPatches);
1376 patchMap_.
setSize(nNewPatches);
1380 fvMeshSubsetPtr_().addFvPatches(newBoundary, syncPar);
1390 const label patchID,
1398 region[iter.key()] = 1;
1407 const label currentRegion,
1408 const bool syncCouples
1412 labelList cellsToRemove(getCellsToRemove(region, currentRegion));
1421 const label currentRegion,
1424 const bool syncCouples
1428 labelList cellsToRemove(getCellsToRemove(region, currentRegion));
1459 pointMap_ = map().pointMap();
1460 faceMap_ = map().faceMap();
1461 cellMap_ = map().cellMap();
1468 return fvMeshSubsetPtr_.valid();
1476 return fvMeshSubsetPtr_();
1484 return fvMeshSubsetPtr_();
1506 if (!faceFlipMapPtr_.valid())
1520 for (
label subFaceI = 0; subFaceI < subInt; subFaceI++)
1522 faceFlipMap[subFaceI] = subToBaseFace[subFaceI]+1;
1524 for (
label subFaceI = subInt; subFaceI < subOwn.
size(); subFaceI++)
1526 label faceI = subToBaseFace[subFaceI];
1527 if (subToBaseCell[subOwn[subFaceI]] == own[faceI])
1529 faceFlipMap[subFaceI] = faceI+1;
1533 faceFlipMap[subFaceI] = -faceI-1;
1538 return faceFlipMapPtr_();
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
static void scatterList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Scatter data. Reverse of gatherList.
const labelList & patchMap() const
Return patch map.
autoPtr< IOobject > clone() const
Clone.
ZoneMesh< faceZone, polyMesh > faceZoneMesh
A ZoneMesh with the type faceZone.
static void listCombineScatter(const List< commsStruct > &comms, List< T > &Value, const int tag, const label comm)
Scatter data. Reverse of combineGather.
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
A face is a list of labels corresponding to mesh vertices.
void setRefinement(const labelList &cellsToRemove, const labelList &facesToExpose, const labelList &patchIDs, polyTopoChange &) const
Play commands into polyTopoChange to remove cells.
autoPtr< mapPolyMesh > makeMesh(autoPtr< fvMesh > &newMesh, const IOobject &io, const polyMesh &mesh, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Create new mesh with old mesh patches.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
label nInternalFaces() const
virtual const labelList & faceNeighbour() const
Return face neighbour.
const labelList & faceFlipMap() const
Return face map with sign to encode flipped faces.
void size(const label)
Override size to be inconsistent with allocated storage.
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Given list of cells to remove insert all the topology changes.
Addressing for all faces on surface of mesh. Can either be read from polyMesh or from triSurface...
Ostream & endl(Ostream &os)
Add newline and flush stream.
label findPatchID(const word &patchName) const
Find patch index given a name.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
const cellList & cells() const
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
label size() const
Return number of elements in table.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
virtual const pointField & points() const
Return raw points.
List< bool > boolList
Bool container classes.
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
labelList getExposedFaces(const labelList ®ion, const label currentRegion, const bool syncCouples=true) const
Two step subsetting.
bool found(const Key &) const
Return true if hashedEntry is found in table.
Xfer< T > xferMove(T &)
Construct by transferring the contents of the arg.
wordList names() const
Return a list of patch names.
wordList patchNames(nPatches)
virtual const labelList & faceOwner() const
Return face owner.
const labelList & faceMap() const
Return face map.
static const label labelMax
List< label > labelList
A List of labels.
ZoneMesh< pointZone, polyMesh > pointZoneMesh
A ZoneMesh with the type pointZone.
virtual const faceList & faces() const
Return raw faces.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
errorManip< error > abort(error &err)
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
Empty front and back plane patch. Used for 2-D geometries.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
word name(const complex &)
Return a string representation of a complex.
label size() const
Return the number of elements in the UPtrList.
void setSize(const label)
Reset size of List.
Template functions to aid in the implementation of demand driven data.
static bool & parRun()
Is this a parallel run?
static label nProcs(const label communicator=0)
Number of processes in parallel run.
const fvMesh & subMesh() const
Return reference to subset mesh.
A cell is defined as a list of faces with extra functionality.
prefixOSstream Pout(cout, "Pout")
Mesh data needed to do the Finite Volume discretisation.
Direct mesh changes based on v1.3 polyTopoChange syntax.
label start() const
Return start label of this patch in the polyMesh face list.
ListType subset(const UList< T > &select, const T &value, const ListType &)
Extract elements of List when select is a certain value.
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
labelList getExposedFaces(const labelList &cellsToRemove) const
Get labels of exposed faces.
List< Key > toc() const
Return the table of contents.
A patch is a list of labels that address the faces in the global face list.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
void setLargeCellSubset(const labelList ®ion, const label currentRegion, const label patchID=-1, const bool syncCouples=true)
Set the subset from all cells with region == currentRegion.
const labelList & cellMap() const
Return cell map.
const fvMesh & baseMesh() const
Original mesh.
ZoneMesh< cellZone, polyMesh > cellZoneMesh
A ZoneMesh with the type cellZone.
void setCellSubset(const labelHashSet &globalCellMap, const label patchID=-1)
Set the subset. Create "oldInternalFaces" patch for exposed.
const labelList & pointMap() const
Return point map.
bool hasSubMesh() const
Have subMesh?