53 void Foam::removeFaces::changeCellRegion
56 const label oldRegion,
57 const label newRegion,
61 if (cellRegion[cellI] == oldRegion)
63 cellRegion[cellI] = newRegion;
71 changeCellRegion(cCells[i], oldRegion, newRegion, cellRegion);
78 Foam::label Foam::removeFaces::changeFaceRegion
84 const label newRegion,
91 if (faceRegion[faceI] == -1 && !removedFace[faceI])
93 faceRegion[faceI] = newRegion;
98 DynamicList<label> fe;
99 DynamicList<label> ef;
104 label edgeI = fEdges[i];
106 if (nFacesPerEdge[edgeI] >= 0 && nFacesPerEdge[edgeI] <= 2)
108 const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
112 label nbrFaceI = eFaces[j];
114 const labelList& fEdges1 = mesh_.faceEdges(nbrFaceI, fe);
116 nChanged += changeFaceRegion
148 boolList affectedFace(mesh_.nFaces(),
false);
153 label region = cellRegion[cellI];
155 if (region != -1 && (cellI != cellRegionMaster[region]))
157 const labelList& cFaces = mesh_.cells()[cellI];
161 affectedFace[cFaces[cFaceI]] =
true;
169 affectedFace[facesToRemove[i]] =
true;
175 const labelList& eFaces = mesh_.edgeFaces(iter.key());
179 affectedFace[eFaces[eFaceI]] =
true;
186 label pointI = iter.key();
192 affectedFace[pFaces[pFaceI]] =
true;
202 const fileName& fName
206 Pout<<
"removeFaces::writeOBJ : Writing faces to file "
209 const pointField& localPoints = fp.localPoints();
216 const faceList& localFaces = fp.localFaces();
220 const face&
f = localFaces[i];
226 str<<
' ' << f[fp]+1;
234 void Foam::removeFaces::mergeFaces
240 polyTopoChange& meshMod
257 if (fp.edgeLoops().size() != 1)
259 writeOBJ(fp, mesh_.time().path()/
"facesToBeMerged.obj");
261 <<
"Cannot merge faces " << faceLabels
262 <<
" into single face since outside vertices " << fp.edgeLoops()
263 <<
" do not form single loop but form " << fp.edgeLoops().size()
267 const labelList& edgeLoop = fp.edgeLoops()[0];
274 label masterIndex = -1;
275 bool reverseLoop =
false;
282 label faceI = pFaces[i];
284 const face&
f = fp.localFaces()[faceI];
286 label index1 =
findIndex(f, edgeLoop[1]);
291 label index0 =
findIndex(f, edgeLoop[0]);
295 if (index1 == f.fcIndex(index0))
301 else if (index1 == f.rcIndex(index0))
311 if (masterIndex == -1)
313 writeOBJ(fp, mesh_.time().path()/
"facesToBeMerged.obj");
323 label faceI = faceLabels[masterIndex];
325 label own = mesh_.faceOwner()[faceI];
327 if (cellRegion[own] != -1)
329 own = cellRegionMaster[cellRegion[own]];
332 label patchID, zoneID, zoneFlip;
334 getFaceInfo(faceI, patchID, zoneID, zoneFlip);
338 if (mesh_.isInternalFace(faceI))
340 nei = mesh_.faceNeighbour()[faceI];
342 if (cellRegion[nei] != -1)
344 nei = cellRegionMaster[cellRegion[nei]];
349 DynamicList<label> faceVerts(edgeLoop.size());
353 label pointI = fp.meshPoints()[edgeLoop[i]];
355 if (pointsToRemove.found(pointI))
362 faceVerts.append(pointI);
367 mergedFace.transfer(faceVerts);
402 forAll(faceLabels, patchFaceI)
404 if (patchFaceI != masterIndex)
408 meshMod.setAction(polyRemoveFace(faceLabels[patchFaceI], faceI));
415 void Foam::removeFaces::getFaceInfo
426 if (!mesh_.isInternalFace(faceI))
428 patchID = mesh_.boundaryMesh().whichPatch(faceI);
431 zoneID = mesh_.faceZones().whichZone(faceI);
437 const faceZone& fZone = mesh_.faceZones()[zoneID];
439 zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
451 const face& f = mesh_.faces()[faceI];
461 if (!pointsToRemove.found(vertI))
463 newFace[newFp++] = vertI;
467 newFace.setSize(newFp);
469 return face(newFace);
474 void Foam::removeFaces::modFace
477 const label masterFaceID,
480 const bool flipFaceFlux,
481 const label newPatchID,
482 const bool removeFromZone,
486 polyTopoChange& meshMod
489 if ((nei == -1) || (own < nei))
561 Foam::removeFaces::removeFaces
588 const labelList& faceOwner = mesh_.faceOwner();
589 const labelList& faceNeighbour = mesh_.faceNeighbour();
591 cellRegion.
setSize(mesh_.nCells());
594 regionMaster.
setSize(mesh_.nCells());
601 label faceI = facesToRemove[i];
603 if (!mesh_.isInternalFace(faceI))
607 "removeFaces::compatibleRemoves(const labelList&"
608 ", labelList&, labelList&, labelList&)"
613 label own = faceOwner[faceI];
614 label nei = faceNeighbour[faceI];
616 label region0 = cellRegion[own];
617 label region1 = cellRegion[nei];
624 cellRegion[own] = nRegions;
625 cellRegion[nei] = nRegions;
628 regionMaster[nRegions] = own;
634 cellRegion[own] = region1;
636 regionMaster[region1] =
min(own, regionMaster[region1]);
644 cellRegion[nei] = region0;
648 else if (region0 != region1)
651 label freedRegion = -1;
652 label keptRegion = -1;
654 if (region0 < region1)
664 keptRegion = region0;
665 freedRegion = region1;
667 else if (region1 < region0)
677 keptRegion = region1;
678 freedRegion = region0;
681 label master0 = regionMaster[region0];
682 label master1 = regionMaster[region1];
684 regionMaster[freedRegion] = -1;
685 regionMaster[keptRegion] =
min(master0, master1);
690 regionMaster.
setSize(nRegions);
701 label r = cellRegion[cellI];
707 if (cellI < regionMaster[r])
711 "removeFaces::compatibleRemoves(const labelList&"
712 ", labelList&, labelList&, labelList&)"
713 ) <<
"Not lowest numbered : cell:" << cellI
715 <<
" regionmaster:" << regionMaster[r]
723 if (nCells[region] == 1)
727 "removeFaces::compatibleRemoves(const labelList&"
728 ", labelList&, labelList&, labelList&)"
729 ) <<
"Region " << region
730 <<
" has only " << nCells[region] <<
" cells in it"
738 label nUsedRegions = 0;
742 if (regionMaster[i] != -1)
751 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
753 label own = faceOwner[faceI];
754 label nei = faceNeighbour[faceI];
756 if (cellRegion[own] != -1 && cellRegion[own] == cellRegion[nei])
760 allFacesToRemove.
append(faceI);
764 newFacesToRemove.
transfer(allFacesToRemove);
780 faceSet facesToRemove(mesh_,
"facesToRemove", faceLabels);
781 Pout<<
"Writing faces to remove to faceSet " << facesToRemove.
name()
783 facesToRemove.
write();
787 boolList removedFace(mesh_.nFaces(),
false);
791 label faceI = faceLabels[i];
793 if (!mesh_.isInternalFace(faceI))
797 "removeFaces::setRefinement(const labelList&"
798 ", const labelList&, const labelList&, polyTopoChange&)"
799 ) <<
"Face to remove is not internal face:" << faceI
803 removedFace[faceI] =
true;
816 labelList faceRegion(mesh_.nFaces(), -1);
831 labelList nFacesPerEdge(mesh_.nEdges(), -1);
836 label faceI = faceLabels[i];
838 const labelList& fEdges = mesh_.faceEdges(faceI, fe);
842 label edgeI = fEdges[i];
844 if (nFacesPerEdge[edgeI] == -1)
846 nFacesPerEdge[edgeI] = mesh_.edgeFaces(edgeI, ef).
size()-1;
850 nFacesPerEdge[edgeI]--;
860 forAll(mesh_.edges(), edgeI)
862 if (nFacesPerEdge[edgeI] == -1)
867 const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
869 if (eFaces.
size() > 2)
871 nFacesPerEdge[edgeI] = eFaces.
size();
873 else if (eFaces.
size() == 2)
879 const edge&
e = mesh_.edges()[edgeI];
882 <<
"Problem : edge has too few face neighbours:"
886 <<
" coords:" << mesh_.points()[e[0]]
887 << mesh_.points()[e[1]]
897 OFstream str(mesh_.time().path()/
"edgesWithTwoFaces.obj");
898 Pout<<
"Dumping edgesWithTwoFaces to " << str.
name() <<
endl;
901 forAll(nFacesPerEdge, edgeI)
903 if (nFacesPerEdge[edgeI] == 2)
906 const edge&
e = mesh_.edges()[edgeI];
912 str<<
"l " << vertI-1 <<
' ' << vertI <<
nl;
923 forAll(nFacesPerEdge, edgeI)
925 if (nFacesPerEdge[edgeI] == 2)
931 const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
935 label faceI = eFaces[i];
937 if (!removedFace[faceI] && !mesh_.isInternalFace(faceI))
951 if (f0 != -1 && f1 != -1)
959 if (patch0 != patch1)
963 <<
"not merging faces " << f0 <<
" and "
964 << f1 <<
" across patch boundary edge " << edgeI
968 nFacesPerEdge[edgeI] = 3;
970 else if (minCos_ < 1 && minCos_ > -1)
980 & n0[f1 - pp0.
start()]
986 <<
"not merging faces " << f0 <<
" and "
987 << f1 <<
" across edge " << edgeI
992 nFacesPerEdge[edgeI] = 3;
996 else if (f0 != -1 || f1 != -1)
998 const edge&
e = mesh_.edges()[edgeI];
1002 <<
"Problem : edge would have one boundary face"
1003 <<
" and one internal face using it." <<
endl
1004 <<
"Your remove pattern is probably incorrect." <<
endl
1006 <<
" nFaces:" << nFacesPerEdge[edgeI]
1007 <<
" vertices:" << e
1008 <<
" coords:" << mesh_.points()[e[0]]
1009 << mesh_.points()[e[1]]
1020 forAll(nFacesPerEdge, edgeI)
1022 if (nFacesPerEdge[edgeI] == 1)
1024 const edge&
e = mesh_.edges()[edgeI];
1027 <<
"Problem : edge would get 1 face using it only"
1028 <<
" edge:" << edgeI
1029 <<
" nFaces:" << nFacesPerEdge[edgeI]
1030 <<
" vertices:" << e
1031 <<
" coords:" << mesh_.points()[e[0]]
1032 <<
' ' << mesh_.points()[e[1]]
1092 forAll(nFacesPerEdge, edgeI)
1094 if (nFacesPerEdge[edgeI] == 0)
1097 edgesToRemove.insert(edgeI);
1099 else if (nFacesPerEdge[edgeI] == 1)
1103 else if (nFacesPerEdge[edgeI] == 2)
1106 edgesToRemove.insert(edgeI);
1112 OFstream str(mesh_.time().path()/
"edgesToRemove.obj");
1113 Pout<<
"Dumping edgesToRemove to " << str.
name() <<
endl;
1119 const edge&
e = mesh_.edges()[iter.key()];
1125 str<<
"l " << vertI-1 <<
' ' << vertI <<
nl;
1133 label startFaceI = 0;
1138 for (; startFaceI < mesh_.nFaces(); startFaceI++)
1140 if (faceRegion[startFaceI] == -1 && !removedFace[startFaceI])
1146 if (startFaceI == mesh_.nFaces())
1154 label nRegion = changeFaceRegion
1161 mesh_.faceEdges(startFaceI, fe),
1169 else if (nRegion == 1)
1172 faceRegion[startFaceI] = -2;
1198 label faceI = mesh_.nInternalFaces();
1199 faceI < mesh_.nFaces();
1204 label nbrRegion = nbrFaceRegion[faceI];
1205 label myRegion = faceRegion[faceI];
1207 if (myRegion <= -1 || nbrRegion <= -1)
1209 if (nbrRegion != myRegion)
1212 <<
"Inconsistent face region across coupled patches."
1214 <<
"This side has for faceI:" << faceI
1215 <<
" region:" << myRegion <<
endl
1216 <<
"The other side has region:" << nbrRegion
1218 <<
"(region -1 means face is to be deleted)"
1222 else if (toNbrRegion[myRegion] == -1)
1225 toNbrRegion[myRegion] = nbrRegion;
1230 if (toNbrRegion[myRegion] != nbrRegion)
1233 <<
"Inconsistent face region across coupled patches."
1235 <<
"This side has for faceI:" << faceI
1236 <<
" region:" << myRegion
1237 <<
" with coupled neighbouring regions:"
1238 << toNbrRegion[myRegion] <<
" and "
1268 labelList nEdgesPerPoint(mesh_.nPoints());
1272 forAll(pointEdges, pointI)
1274 nEdgesPerPoint[pointI] = pointEdges[pointI].size();
1280 const edge&
e = mesh_.edges()[iter.key()];
1284 nEdgesPerPoint[e[i]]--;
1289 forAll(nEdgesPerPoint, pointI)
1291 if (nEdgesPerPoint[pointI] == 1)
1294 <<
"Problem : point would get 1 edge using it only."
1295 <<
" pointI:" << pointI
1296 <<
" coord:" << mesh_.points()[pointI]
1312 forAll(nEdgesPerPoint, pointI)
1314 if (nEdgesPerPoint[pointI] == 0)
1316 pointsToRemove.insert(pointI);
1318 else if (nEdgesPerPoint[pointI] == 1)
1322 else if (nEdgesPerPoint[pointI] == 2)
1325 pointsToRemove.insert(pointI);
1333 OFstream str(mesh_.time().path()/
"pointsToRemove.obj");
1334 Pout<<
"Dumping pointsToRemove to " << str.
name() <<
endl;
1378 label faceI = faceLabels[
labelI];
1382 if (affectedFace[faceI])
1384 affectedFace[faceI] =
false;
1394 label pointI = iter.key();
1401 forAll(cellRegion, cellI)
1403 label region = cellRegion[cellI];
1405 if (region != -1 && (cellI != cellRegionMaster[region]))
1420 forAll(regionToFaces, regionI)
1422 const labelList& rFaces = regionToFaces[regionI];
1424 if (rFaces.
size() <= 1)
1427 <<
"Region:" << regionI
1428 <<
" contains only faces " << rFaces
1444 affectedFace[rFaces[i]] =
false;
1455 forAll(affectedFace, faceI)
1457 if (affectedFace[faceI])
1459 affectedFace[faceI] =
false;
1461 face f(filterFace(pointsToRemove, faceI));
1463 label own = mesh_.faceOwner()[faceI];
1465 if (cellRegion[own] != -1)
1467 own = cellRegionMaster[cellRegion[own]];
1470 label patchID, zoneID, zoneFlip;
1472 getFaceInfo(faceI, patchID, zoneID, zoneFlip);
1476 if (mesh_.isInternalFace(faceI))
1478 nei = mesh_.faceNeighbour()[faceI];
1480 if (cellRegion[nei] != -1)
1482 nei = cellRegionMaster[cellRegion[nei]];