48 bool Foam::fvMeshSubset::checkCellSubset()
const
50 if (fvMeshSubsetPtr_.empty())
52 FatalErrorIn(
"bool fvMeshSubset::checkCellSubset() const")
53 <<
"Mesh subset not set. Please set the cell map using "
54 <<
"void setCellSubset(const labelHashSet& cellsToSubset)" <<
endl
55 <<
"before attempting to access subset data"
67 void Foam::fvMeshSubset::markPoints
76 pointMap.insert(curPoints[pointI], 0);
81 void Foam::fvMeshSubset::markPoints
89 pointMap[curPoints[pointI]] = 0;
96 void Foam::fvMeshSubset::doCoupledPatches
102 const polyBoundaryMesh& oldPatches = baseMesh().boundaryMesh();
104 label nUncoupled = 0;
109 forAll (oldPatches, oldPatchI)
111 const polyPatch& pp = oldPatches[oldPatchI];
113 if (isA<processorPolyPatch>(pp))
115 const processorPolyPatch& procPatch =
116 refCast<const processorPolyPatch>(pp);
121 procPatch.neighbProcNo()
125 << SubList<label>(nCellsUsingFace, pp.size(), pp.start());
131 forAll (oldPatches, oldPatchI)
133 const polyPatch& pp = oldPatches[oldPatchI];
135 if (isA<processorPolyPatch>(pp))
137 const processorPolyPatch& procPatch =
138 refCast<const processorPolyPatch>(pp);
140 IPstream fromNeighbour
143 procPatch.neighbProcNo()
146 labelList nbrCellsUsingFace(fromNeighbour);
154 nCellsUsingFace[pp.start()+i] == 1
155 && nbrCellsUsingFace[i] == 0
160 nCellsUsingFace[pp.start()+i] = 3;
169 forAll (oldPatches, oldPatchI)
171 const polyPatch& pp = oldPatches[oldPatchI];
173 if (isA<cyclicPolyPatch>(pp))
175 const cyclicPolyPatch& cycPatch =
176 refCast<const cyclicPolyPatch>(pp);
180 label thisFaceI = cycPatch.start() + i;
181 label otherFaceI = cycPatch.transformGlobalFace(thisFaceI);
185 nCellsUsingFace[thisFaceI] == 1
186 && nCellsUsingFace[otherFaceI] == 0
189 nCellsUsingFace[thisFaceI] = 3;
198 reduce(nUncoupled, sumOp<label>());
203 Info<<
"Uncoupled " << nUncoupled <<
" faces on coupled patches. "
204 <<
"(processorPolyPatch, cyclicPolyPatch)" <<
nl
205 <<
"You might need to run couplePatches to restore the patch face"
206 <<
" ordering." <<
nl <<
endl;
220 forAll(selectedElements, i)
222 selected[selectedElements[i]] =
true;
229 if (selected[subsetMap[i]])
241 if (selected[subsetMap[i]])
243 subsettedElements[n++] = i;
247 return subsettedElements;
251 void Foam::fvMeshSubset::subsetZones()
258 List<pointZone*> pZonePtrs(pointZones.size());
262 const pointZone& pz = pointZones[i];
264 pZonePtrs[i] =
new pointZone
269 fvMeshSubsetPtr_().pointZones()
281 List<faceZone*> fZonePtrs(faceZones.size());
285 const faceZone& fz = faceZones[i];
308 if (zone[faceMap()[j]] != 0)
316 forAll(faceMap(), subFaceI)
318 label meshFaceI = faceMap()[subFaceI];
319 if (zone[meshFaceI] != 0)
321 subAddressing[nSub] = subFaceI;
322 label subOwner = subMesh().faceOwner()[subFaceI];
323 label baseOwner = baseMesh().faceOwner()[meshFaceI];
325 bool sameOwner = (cellMap()[subOwner] == baseOwner);
326 bool flip = (zone[meshFaceI] == 1);
327 subFlipStatus[nSub] = (sameOwner == flip);
333 fZonePtrs[i] =
new faceZone
339 fvMeshSubsetPtr_().faceZones()
346 List<cellZone*> cZonePtrs(cellZones.size());
350 const cellZone& cz = cellZones[i];
352 cZonePtrs[i] =
new cellZone
355 subset(baseMesh().nCells(), cz, cellMap()),
357 fvMeshSubsetPtr_().cellZones()
363 fvMeshSubsetPtr_().addZones(pZonePtrs, fZonePtrs, cZonePtrs);
370 Foam::fvMeshSubset::fvMeshSubset(
const fvMesh& baseMesh)
373 fvMeshSubsetPtr_(NULL),
391 const cellList& oldCells = baseMesh().cells();
392 const faceList& oldFaces = baseMesh().faces();
393 const pointField& oldPoints = baseMesh().points();
394 const labelList& oldOwner = baseMesh().faceOwner();
395 const labelList& oldNeighbour = baseMesh().faceNeighbour();
397 label wantedPatchID = patchID;
399 if (wantedPatchID == -1)
403 wantedPatchID = oldPatches.
findPatchID(
"oldInternalFaces");
405 else if (wantedPatchID < 0 || wantedPatchID >= oldPatches.
size())
409 "fvMeshSubset::setCellSubset(const labelHashSet&"
410 ", const label patchID)"
411 ) <<
"Non-existing patch index " << wantedPatchID <<
endl
412 <<
"Should be between 0 and " << oldPatches.
size()-1
417 cellMap_ = globalCellMap.
toc();
423 const label avgNFacesPerCell = 6;
424 const label avgNPointsPerFace = 4;
427 label nCellsInSet = cellMap_.size();
431 Map<label> facesToSubset(avgNFacesPerCell*nCellsInSet);
436 const labelList& curFaces = oldCells[cellMap_[cellI]];
440 if (!facesToSubset.
found(curFaces[faceI]))
442 facesToSubset.
insert(curFaces[faceI], 1);
446 facesToSubset[curFaces[faceI]]++;
455 Map<label> globalPointMap(avgNPointsPerFace*facesToSubset.
size());
462 faceMap_.setSize(facesToc.
size());
467 if (facesToSubset[facesToc[faceI]] == 2)
470 faceMap_[globalFaceMap.size()] = facesToc[faceI];
471 globalFaceMap.
insert(facesToc[faceI], globalFaceMap.
size());
474 markPoints(oldFaces[facesToc[faceI]], globalPointMap);
479 label nInternalFaces = globalFaceMap.size();
483 label oldPatchStart = labelMax;
484 if (wantedPatchID != -1)
486 oldPatchStart = oldPatches[wantedPatchID].start();
493 for (; faceI< facesToc.
size(); faceI++)
495 if (facesToc[faceI] >= oldPatchStart)
501 !baseMesh().isInternalFace(facesToc[faceI])
502 && facesToSubset[facesToc[faceI]] == 1
506 faceMap_[globalFaceMap.size()] = facesToc[faceI];
507 globalFaceMap.insert(facesToc[faceI], globalFaceMap.size());
510 markPoints(oldFaces[facesToc[faceI]], globalPointMap);
515 forAll(facesToc, intFaceI)
519 baseMesh().isInternalFace(facesToc[intFaceI])
520 && facesToSubset[facesToc[intFaceI]] == 1
524 faceMap_[globalFaceMap.size()] = facesToc[intFaceI];
525 globalFaceMap.
insert(facesToc[intFaceI], globalFaceMap.
size());
528 markPoints(oldFaces[facesToc[intFaceI]], globalPointMap);
533 for (; faceI< facesToc.
size(); faceI++)
537 !baseMesh().isInternalFace(facesToc[faceI])
538 && facesToSubset[facesToc[faceI]] == 1
542 faceMap_[globalFaceMap.size()] = facesToc[faceI];
543 globalFaceMap.insert(facesToc[faceI], globalFaceMap.size());
546 markPoints(oldFaces[facesToc[faceI]], globalPointMap);
553 pointMap_ = globalPointMap.toc();
556 forAll (pointMap_, pointI)
558 globalPointMap[pointMap_[pointI]] = pointI;
561 Pout <<
"Number of cells in new mesh: " << nCellsInSet <<
endl;
562 Pout <<
"Number of faces in new mesh: " << globalFaceMap.size() <<
endl;
563 Pout <<
"Number of points in new mesh: " << globalPointMap.size() <<
endl;
568 label nNewPoints = 0;
570 forAll (pointMap_, pointI)
572 newPoints[nNewPoints] = oldPoints[pointMap_[pointI]];
576 faceList newFaces(globalFaceMap.size());
581 for (label faceI = 0; faceI < nInternalFaces; faceI++)
583 const face& oldF = oldFaces[faceMap_[faceI]];
589 newF[i] = globalPointMap[oldF[i]];
592 newFaces[nNewFaces] = newF;
598 label nbSize = oldPatches.
size();
599 label oldInternalPatchID = -1;
601 if (wantedPatchID == -1)
605 oldInternalPatchID = nbSize;
611 oldInternalPatchID = wantedPatchID;
621 for (label faceI = nInternalFaces; faceI < faceMap_.size(); faceI++)
623 label oldFaceI = faceMap_[faceI];
625 face oldF = oldFaces[oldFaceI];
628 if (baseMesh().isInternalFace(oldFaceI))
633 !globalCellMap.
found(oldOwner[oldFaceI])
634 && globalCellMap.
found(oldNeighbour[oldFaceI])
637 oldF = oldFaces[oldFaceI].reverseFace();
641 boundaryPatchSizes[oldInternalPatchID]++;
646 label patchOfFace = oldPatches.
whichPatch(oldFaceI);
649 boundaryPatchSizes[patchOfFace]++;
656 newF[i] = globalPointMap[oldF[i]];
659 newFaces[nNewFaces] = newF;
672 const labelList& oldC = oldCells[cellMap_[cellI]];
678 newC[i] = globalFaceMap[oldC[i]];
681 newCells[nNewCells] =
cell(newC);
687 fvMeshSubsetPtr_.clear();
689 fvMeshSubsetPtr_.reset
695 baseMesh().
name() +
"SubSet",
710 patchMap_.setSize(nbSize);
711 label nNewPatches = 0;
712 label patchStart = nInternalFaces;
714 forAll(oldPatches, patchI)
716 if (boundaryPatchSizes[patchI] > 0)
719 newBoundary[nNewPatches] = oldPatches[patchI].
clone
723 boundaryPatchSizes[patchI],
727 patchStart += boundaryPatchSizes[patchI];
728 patchMap_[nNewPatches] = patchI;
733 if (wantedPatchID == -1)
736 if (boundaryPatchSizes[oldInternalPatchID] > 0)
741 boundaryPatchSizes[oldInternalPatchID],
749 patchMap_[nNewPatches] = -1;
755 newBoundary.
setSize(nNewPatches);
756 patchMap_.setSize(nNewPatches);
759 fvMeshSubsetPtr_().addFvPatches(newBoundary);
769 const label currentRegion,
774 const cellList& oldCells = baseMesh().cells();
775 const faceList& oldFaces = baseMesh().faces();
776 const pointField& oldPoints = baseMesh().points();
777 const labelList& oldOwner = baseMesh().faceOwner();
778 const labelList& oldNeighbour = baseMesh().faceNeighbour();
780 const label oldNInternalFaces = baseMesh().nInternalFaces();
784 if (region.
size() != oldCells.
size())
788 "fvMeshSubset::setCellSubset(const labelList&"
789 ", const label, const label, const bool)"
790 ) <<
"Size of region " << region.
size()
791 <<
" is not equal to number of cells in mesh " << oldCells.
size()
796 label wantedPatchID = patchID;
798 if (wantedPatchID == -1)
802 wantedPatchID = oldPatches.
findPatchID(
"oldInternalFaces");
804 else if (wantedPatchID < 0 || wantedPatchID >= oldPatches.
size())
808 "fvMeshSubset::setCellSubset(const labelList&"
809 ", const label, const label, const bool)"
810 ) <<
"Non-existing patch index " << wantedPatchID <<
endl
811 <<
"Should be between 0 and " << oldPatches.
size()-1
817 cellMap_.setSize(oldCells.
size());
818 label nCellsInSet = 0;
822 if (region[oldCellI] == currentRegion)
824 cellMap_[nCellsInSet++] = oldCellI;
827 cellMap_.setSize(nCellsInSet);
842 label nFacesInSet = 0;
843 forAll (oldFaces, oldFaceI)
845 bool faceUsed =
false;
847 if (region[oldOwner[oldFaceI]] == currentRegion)
849 nCellsUsingFace[oldFaceI]++;
855 baseMesh().isInternalFace(oldFaceI)
856 && (region[oldNeighbour[oldFaceI]] == currentRegion)
859 nCellsUsingFace[oldFaceI]++;
868 faceMap_.setSize(nFacesInSet);
871 doCoupledPatches(syncPar, nCellsUsingFace);
875 label oldInternalPatchID = 0;
878 label nextPatchID = oldPatches.
size();
884 label nbSize = oldPatches.
size();
886 if (wantedPatchID == -1)
892 forAll(oldPatches, patchI)
894 if (isA<processorPolyPatch>(oldPatches[patchI]))
896 nextPatchID = patchI;
899 oldInternalPatchID++;
905 for (label oldPatchI = 0; oldPatchI < nextPatchID; oldPatchI++)
907 globalPatchMap[oldPatchI] = oldPatchI;
911 label oldPatchI = nextPatchID;
912 oldPatchI < oldPatches.
size();
916 globalPatchMap[oldPatchI] = oldPatchI+1;
921 oldInternalPatchID = wantedPatchID;
922 nextPatchID = wantedPatchID+1;
938 for (label oldFaceI = 0; oldFaceI < oldNInternalFaces; oldFaceI++)
940 if (nCellsUsingFace[oldFaceI] == 2)
942 globalFaceMap[oldFaceI] = faceI;
943 faceMap_[faceI++] = oldFaceI;
946 markPoints(oldFaces[oldFaceI], globalPointMap);
951 label nInternalFaces = faceI;
957 oldPatchI < oldPatches.
size()
958 && oldPatchI < nextPatchID;
962 const polyPatch& oldPatch = oldPatches[oldPatchI];
964 label oldFaceI = oldPatch.
start();
968 if (nCellsUsingFace[oldFaceI] == 1)
973 globalFaceMap[oldFaceI] = faceI;
974 faceMap_[faceI++] = oldFaceI;
977 markPoints(oldFaces[oldFaceI], globalPointMap);
980 boundaryPatchSizes[globalPatchMap[oldPatchI]]++;
987 for (label oldFaceI = 0; oldFaceI < oldNInternalFaces; oldFaceI++)
989 if (nCellsUsingFace[oldFaceI] == 1)
991 globalFaceMap[oldFaceI] = faceI;
992 faceMap_[faceI++] = oldFaceI;
995 markPoints(oldFaces[oldFaceI], globalPointMap);
998 boundaryPatchSizes[oldInternalPatchID]++;
1005 label oldFaceI = oldNInternalFaces;
1006 oldFaceI < oldFaces.
size();
1010 if (nCellsUsingFace[oldFaceI] == 3)
1012 globalFaceMap[oldFaceI] = faceI;
1013 faceMap_[faceI++] = oldFaceI;
1016 markPoints(oldFaces[oldFaceI], globalPointMap);
1019 boundaryPatchSizes[oldInternalPatchID]++;
1026 label oldPatchI = nextPatchID;
1027 oldPatchI < oldPatches.
size();
1031 const polyPatch& oldPatch = oldPatches[oldPatchI];
1033 label oldFaceI = oldPatch.
start();
1037 if (nCellsUsingFace[oldFaceI] == 1)
1042 globalFaceMap[oldFaceI] = faceI;
1043 faceMap_[faceI++] = oldFaceI;
1046 markPoints(oldFaces[oldFaceI], globalPointMap);
1049 boundaryPatchSizes[globalPatchMap[oldPatchI]]++;
1055 if (faceI != nFacesInSet)
1059 "fvMeshSubset::setCellSubset(const labelList&"
1060 ", const label, const label, const bool)"
1066 label nPointsInSet = 0;
1068 forAll(globalPointMap, pointI)
1070 if (globalPointMap[pointI] != -1)
1075 pointMap_.setSize(nPointsInSet);
1079 forAll(globalPointMap, pointI)
1081 if (globalPointMap[pointI] != -1)
1083 pointMap_[nPointsInSet] = pointI;
1084 globalPointMap[pointI] = nPointsInSet;
1096 label nNewPoints = 0;
1098 forAll (pointMap_, pointI)
1100 newPoints[nNewPoints] = oldPoints[pointMap_[pointI]];
1104 faceList newFaces(faceMap_.size());
1106 label nNewFaces = 0;
1109 for (label faceI = 0; faceI < nInternalFaces; faceI++)
1111 const face& oldF = oldFaces[faceMap_[faceI]];
1117 newF[i] = globalPointMap[oldF[i]];
1120 newFaces[nNewFaces] = newF;
1127 for (label faceI = nInternalFaces; faceI < faceMap_.size(); faceI++)
1129 label oldFaceI = faceMap_[faceI];
1131 face oldF = oldFaces[oldFaceI];
1134 if (baseMesh().isInternalFace(oldFaceI))
1139 region[oldOwner[oldFaceI]] != currentRegion
1140 && region[oldNeighbour[oldFaceI]] == currentRegion
1143 oldF = oldFaces[oldFaceI].reverseFace();
1152 newF[i] = globalPointMap[oldF[i]];
1155 newFaces[nNewFaces] = newF;
1164 label nNewCells = 0;
1168 const labelList& oldC = oldCells[cellMap_[cellI]];
1174 newC[i] = globalFaceMap[oldC[i]];
1177 newCells[nNewCells] =
cell(newC);
1187 fvMeshSubsetPtr_.reset
1208 patchMap_.setSize(nbSize);
1209 label nNewPatches = 0;
1210 label patchStart = nInternalFaces;
1218 labelList globalPatchSizes(boundaryPatchSizes);
1219 globalPatchSizes.
setSize(nextPatchID);
1239 bool samePatches =
true;
1241 for (label procI = 1; procI < patchNames.
size(); procI++)
1243 if (patchNames[procI] != patchNames[0])
1245 samePatches =
false;
1254 globalPatchSizes = labelMax;
1263 label oldPatchI = 0;
1264 oldPatchI < oldPatches.
size()
1265 && oldPatchI < nextPatchID;
1269 label newSize = boundaryPatchSizes[globalPatchMap[oldPatchI]];
1272 newBoundary[nNewPatches] = oldPatches[oldPatchI].
clone
1280 patchStart += newSize;
1281 patchMap_[nNewPatches] = oldPatchI;
1287 if (wantedPatchID == -1)
1289 label oldInternalSize = boundaryPatchSizes[oldInternalPatchID];
1297 if (oldInternalSize > 0)
1302 boundaryPatchSizes[oldInternalPatchID],
1313 patchStart += boundaryPatchSizes[oldInternalPatchID];
1314 patchMap_[nNewPatches] = -1;
1323 label oldPatchI = nextPatchID;
1324 oldPatchI < oldPatches.
size();
1328 label newSize = boundaryPatchSizes[globalPatchMap[oldPatchI]];
1331 newBoundary[nNewPatches] = oldPatches[oldPatchI].
clone
1342 patchStart += newSize;
1343 patchMap_[nNewPatches] = oldPatchI;
1349 newBoundary.
setSize(nNewPatches);
1350 patchMap_.setSize(nNewPatches);
1354 fvMeshSubsetPtr_().addFvPatches(newBoundary, syncPar);
1364 const label patchID,
1368 labelList region(baseMesh().nCells(), 0);
1372 region[iter.key()] = 1;
1374 setLargeCellSubset(region, 1, patchID, syncPar);
1382 return fvMeshSubsetPtr_();
1390 return fvMeshSubsetPtr_();