53 void Foam::removeFaces::changeCellRegion
56 const label oldRegion,
57 const label newRegion,
61 if (cellRegion[celli] == oldRegion)
63 cellRegion[celli] = newRegion;
67 const labelList& cCells = mesh_.cellCells()[celli];
71 changeCellRegion(cCells[i], oldRegion, newRegion, cellRegion);
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();
188 const labelList& pFaces = mesh_.pointFaces()[pointi];
192 affectedFace[pFaces[pFacei]] =
true;
199 void Foam::removeFaces::writeOBJ
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;
277 const labelList& pFaces = fp.pointFaces()[edgeLoop[0]];
282 label facei = pFaces[i];
284 const face& f = fp.localFaces()[facei];
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))
610 label own = faceOwner[facei];
611 label nei = faceNeighbour[facei];
613 label region0 = cellRegion[own];
614 label region1 = cellRegion[nei];
621 cellRegion[own] = nRegions;
622 cellRegion[nei] = nRegions;
625 regionMaster[nRegions] = own;
631 cellRegion[own] = region1;
633 regionMaster[region1] =
min(own, regionMaster[region1]);
641 cellRegion[nei] = region0;
645 else if (region0 != region1)
648 label freedRegion = -1;
649 label keptRegion = -1;
651 if (region0 < region1)
661 keptRegion = region0;
662 freedRegion = region1;
664 else if (region1 < region0)
674 keptRegion = region1;
675 freedRegion = region0;
678 label master0 = regionMaster[region0];
679 label master1 = regionMaster[region1];
681 regionMaster[freedRegion] = -1;
682 regionMaster[keptRegion] =
min(master0, master1);
687 regionMaster.
setSize(nRegions);
698 label r = cellRegion[celli];
704 if (celli < regionMaster[r])
707 <<
"Not lowest numbered : cell:" << celli
709 <<
" regionmaster:" << regionMaster[r]
717 if (nCells[region] == 1)
720 <<
"Region " << region
721 <<
" has only " << nCells[region] <<
" cells in it" 729 label nUsedRegions = 0;
733 if (regionMaster[i] != -1)
742 for (
label facei = 0; facei < mesh_.nInternalFaces(); facei++)
744 label own = faceOwner[facei];
745 label nei = faceNeighbour[facei];
747 if (cellRegion[own] != -1 && cellRegion[own] == cellRegion[nei])
751 allFacesToRemove.
append(facei);
755 newFacesToRemove.
transfer(allFacesToRemove);
771 faceSet facesToRemove(mesh_,
"facesToRemove", faceLabels);
772 Pout<<
"Writing faces to remove to faceSet " << facesToRemove.
name()
774 facesToRemove.
write();
778 boolList removedFace(mesh_.nFaces(),
false);
782 label facei = faceLabels[i];
784 if (!mesh_.isInternalFace(facei))
787 <<
"Face to remove is not internal face:" << facei
791 removedFace[facei] =
true;
804 labelList faceRegion(mesh_.nFaces(), -1);
819 labelList nFacesPerEdge(mesh_.nEdges(), -1);
824 label facei = faceLabels[i];
826 const labelList& fEdges = mesh_.faceEdges(facei, fe);
830 label edgeI = fEdges[i];
832 if (nFacesPerEdge[edgeI] == -1)
834 nFacesPerEdge[edgeI] = mesh_.edgeFaces(edgeI, ef).
size()-1;
838 nFacesPerEdge[edgeI]--;
848 forAll(mesh_.edges(), edgeI)
850 if (nFacesPerEdge[edgeI] == -1)
855 const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
857 if (eFaces.
size() > 2)
859 nFacesPerEdge[edgeI] = eFaces.
size();
861 else if (eFaces.
size() == 2)
867 const edge&
e = mesh_.edges()[edgeI];
870 <<
"Problem : edge has too few face neighbours:" 874 <<
" coords:" << mesh_.points()[e[0]]
875 << mesh_.points()[e[1]]
885 OFstream str(mesh_.time().path()/
"edgesWithTwoFaces.obj");
886 Pout<<
"Dumping edgesWithTwoFaces to " << str.
name() <<
endl;
889 forAll(nFacesPerEdge, edgeI)
891 if (nFacesPerEdge[edgeI] == 2)
894 const edge&
e = mesh_.edges()[edgeI];
900 str<<
"l " << vertI-1 <<
' ' << vertI <<
nl;
911 forAll(nFacesPerEdge, edgeI)
913 if (nFacesPerEdge[edgeI] == 2)
919 const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
923 label facei = eFaces[i];
925 if (!removedFace[facei] && !mesh_.isInternalFace(facei))
939 if (f0 != -1 && f1 != -1)
947 if (patch0 != patch1)
951 <<
"not merging faces " << f0 <<
" and " 952 << f1 <<
" across patch boundary edge " << edgeI
956 nFacesPerEdge[edgeI] = 3;
958 else if (minCos_ < 1 && minCos_ > -1)
968 & n0[f1 - pp0.
start()]
974 <<
"not merging faces " << f0 <<
" and " 975 << f1 <<
" across edge " << edgeI
980 nFacesPerEdge[edgeI] = 3;
984 else if (f0 != -1 || f1 != -1)
986 const edge&
e = mesh_.edges()[edgeI];
990 <<
"Problem : edge would have one boundary face" 991 <<
" and one internal face using it." <<
endl 992 <<
"Your remove pattern is probably incorrect." <<
endl 994 <<
" nFaces:" << nFacesPerEdge[edgeI]
996 <<
" coords:" << mesh_.points()[e[0]]
997 << mesh_.points()[e[1]]
1008 forAll(nFacesPerEdge, edgeI)
1010 if (nFacesPerEdge[edgeI] == 1)
1012 const edge&
e = mesh_.edges()[edgeI];
1015 <<
"Problem : edge would get 1 face using it only" 1016 <<
" edge:" << edgeI
1017 <<
" nFaces:" << nFacesPerEdge[edgeI]
1018 <<
" vertices:" << e
1019 <<
" coords:" << mesh_.points()[e[0]]
1020 <<
' ' << mesh_.points()[e[1]]
1079 forAll(nFacesPerEdge, edgeI)
1081 if (nFacesPerEdge[edgeI] == 0)
1084 edgesToRemove.insert(edgeI);
1086 else if (nFacesPerEdge[edgeI] == 1)
1090 else if (nFacesPerEdge[edgeI] == 2)
1093 edgesToRemove.insert(edgeI);
1099 OFstream str(mesh_.time().path()/
"edgesToRemove.obj");
1100 Pout<<
"Dumping edgesToRemove to " << str.
name() <<
endl;
1106 const edge&
e = mesh_.edges()[iter.key()];
1112 str<<
"l " << vertI-1 <<
' ' << vertI <<
nl;
1120 label startFacei = 0;
1125 for (; startFacei < mesh_.nFaces(); startFacei++)
1127 if (faceRegion[startFacei] == -1 && !removedFace[startFacei])
1133 if (startFacei == mesh_.nFaces())
1141 label nRegion = changeFaceRegion
1148 mesh_.faceEdges(startFacei, fe),
1156 else if (nRegion == 1)
1159 faceRegion[startFacei] = -2;
1184 label facei = mesh_.nInternalFaces();
1185 facei < mesh_.nFaces();
1190 label nbrRegion = nbrFaceRegion[facei];
1191 label myRegion = faceRegion[facei];
1193 if (myRegion <= -1 || nbrRegion <= -1)
1195 if (nbrRegion != myRegion)
1198 <<
"Inconsistent face region across coupled patches." 1200 <<
"This side has for facei:" << facei
1201 <<
" region:" << myRegion <<
endl 1202 <<
"The other side has region:" << nbrRegion
1204 <<
"(region -1 means face is to be deleted)" 1208 else if (toNbrRegion[myRegion] == -1)
1211 toNbrRegion[myRegion] = nbrRegion;
1216 if (toNbrRegion[myRegion] != nbrRegion)
1219 <<
"Inconsistent face region across coupled patches." 1221 <<
"This side has for facei:" << facei
1222 <<
" region:" << myRegion
1223 <<
" with coupled neighbouring regions:" 1224 << toNbrRegion[myRegion] <<
" and " 1254 labelList nEdgesPerPoint(mesh_.nPoints());
1258 forAll(pointEdges, pointi)
1260 nEdgesPerPoint[pointi] = pointEdges[pointi].size();
1266 const edge&
e = mesh_.edges()[iter.key()];
1270 nEdgesPerPoint[e[i]]--;
1275 forAll(nEdgesPerPoint, pointi)
1277 if (nEdgesPerPoint[pointi] == 1)
1280 <<
"Problem : point would get 1 edge using it only." 1281 <<
" pointi:" << pointi
1282 <<
" coord:" << mesh_.points()[pointi]
1297 forAll(nEdgesPerPoint, pointi)
1299 if (nEdgesPerPoint[pointi] == 0)
1301 pointsToRemove.insert(pointi);
1303 else if (nEdgesPerPoint[pointi] == 1)
1307 else if (nEdgesPerPoint[pointi] == 2)
1310 pointsToRemove.insert(pointi);
1318 OFstream str(mesh_.time().path()/
"pointsToRemove.obj");
1319 Pout<<
"Dumping pointsToRemove to " << str.
name() <<
endl;
1367 if (affectedFace[facei])
1369 affectedFace[facei] =
false;
1379 label pointi = iter.key();
1386 forAll(cellRegion, celli)
1388 label region = cellRegion[celli];
1390 if (region != -1 && (celli != cellRegionMaster[region]))
1405 forAll(regionToFaces, regionI)
1407 const labelList& rFaces = regionToFaces[regionI];
1409 if (rFaces.
size() <= 1)
1412 <<
"Region:" << regionI
1413 <<
" contains only faces " << rFaces
1429 affectedFace[rFaces[i]] =
false;
1440 forAll(affectedFace, facei)
1442 if (affectedFace[facei])
1444 affectedFace[facei] =
false;
1446 face f(filterFace(pointsToRemove, facei));
1448 label own = mesh_.faceOwner()[facei];
1450 if (cellRegion[own] != -1)
1452 own = cellRegionMaster[cellRegion[own]];
1455 label patchID, zoneID, zoneFlip;
1457 getFaceInfo(facei, patchID, zoneID, zoneFlip);
1461 if (mesh_.isInternalFace(facei))
1463 nei = mesh_.faceNeighbour()[facei];
1465 if (cellRegion[nei] != -1)
1467 nei = cellRegionMaster[cellRegion[nei]];
Class containing data for face removal.
#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.
const double e
Elementary charge.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
labelListList invertOneToMany(const label len, const labelUList &)
Invert one-to-many map. Unmapped elements will be size 0.
void size(const label)
Override size to be inconsistent with allocated storage.
Ostream & endl(Ostream &os)
Add newline and flush stream.
PrimitivePatch< face, IndirectList, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
List< bool > boolList
Bool container classes.
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
vectorField pointField
pointField is a vectorField.
Class containing data for cell removal.
static const labelSphericalTensor labelI(1)
Identity labelTensor.
label start() const
Return start label of this patch in the polyMesh face list.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
label compatibleRemoves(const labelList &inPiercedFaces, labelList &cellRegion, labelList &cellRegionMaster, labelList &outPiercedFaces) const
Find faces including those with cells which have the same mastercell.
virtual const fileName & name() const
Return the name of the stream.
void append(const T &)
Append an element at the end of the list.
List< label > labelList
A List of labels.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
errorManip< error > abort(error &err)
prefixOSstream Pout(cout,"Pout")
void reverse(UList< T > &, const label n)
defineTypeNameAndDebug(combustionModel, 0)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void setSize(const label)
Reset size of List.
const Field< PointType > & faceNormals() const
Return face normals for patch.
#define WarningInFunction
Report a warning using Foam::Warning.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
Direct mesh changes based on v1.3 polyTopoChange syntax.
virtual bool write() const
Write using setting from DB.
void setRefinement(const labelList &piercedFaces, const labelList &cellRegion, const labelList &cellRegionMaster, polyTopoChange &) const
Play commands into polyTopoChange to remove faces.
dimensioned< scalar > mag(const dimensioned< Type > &)
Mesh consisting of general polyhedral cells.
A patch is a list of labels that address the faces in the global face list.
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Class containing data for point removal.
static const label labelMin
const word & name() const
Return name.
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.