56 void Foam::removeFaces::changeCellRegion
59 const label oldRegion,
60 const label newRegion,
64 if (cellRegion[celli] == oldRegion)
66 cellRegion[celli] = newRegion;
74 changeCellRegion(cCells[i], oldRegion, newRegion, cellRegion);
81 Foam::label Foam::removeFaces::changeFaceRegion
87 const label newRegion,
94 if (faceRegion[facei] == -1 && !removedFace[facei])
96 faceRegion[facei] = newRegion;
101 DynamicList<label> fe;
102 DynamicList<label> ef;
107 label edgeI = fEdges[i];
109 if (nFacesPerEdge[edgeI] >= 0 && nFacesPerEdge[edgeI] <= 2)
111 const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
115 label nbrFacei = eFaces[j];
117 const labelList& fEdges1 = mesh_.faceEdges(nbrFacei, fe);
119 nChanged += changeFaceRegion
151 boolList affectedFace(mesh_.nFaces(),
false);
156 label region = cellRegion[celli];
158 if (region != -1 && (celli != cellRegionMaster[region]))
160 const labelList& cFaces = mesh_.cells()[celli];
164 affectedFace[cFaces[cFacei]] =
true;
172 affectedFace[facesToRemove[i]] =
true;
176 for (
const label edgei : edgesToRemove)
178 const labelList& eFaces = mesh_.edgeFaces(edgei);
182 affectedFace[eFaces[eFacei]] =
true;
187 for (
const label pointi : pointsToRemove)
193 affectedFace[
pFaces[pFacei]] =
true;
200 void Foam::removeFaces::writeOBJ
203 const fileName& fName
207 Pout<<
"removeFaces::writeOBJ : Writing faces to file " 210 const pointField& localPoints = fp.localPoints();
217 const faceList& localFaces = fp.localFaces();
221 const face&
f = localFaces[i];
227 str<<
' ' <<
f[fp]+1;
235 void Foam::removeFaces::mergeFaces
241 polyTopoChange& meshMod
258 if (fp.edgeLoops().size() != 1)
260 writeOBJ(fp, mesh_.time().path()/
"facesToBeMerged.obj");
263 <<
" into single face since outside vertices " << fp.edgeLoops()
264 <<
" do not form single loop but form " << fp.edgeLoops().
size()
268 const labelList& edgeLoop = fp.edgeLoops()[0];
275 label masterIndex = -1;
276 bool reverseLoop =
false;
285 const face&
f = fp.localFaces()[facei];
287 label index1 =
f.
find(edgeLoop[1]);
292 label index0 =
f.
find(edgeLoop[0]);
302 else if (index1 ==
f.
rcIndex(index0))
312 if (masterIndex == -1)
314 writeOBJ(fp, mesh_.time().path()/
"facesToBeMerged.obj");
326 label own = mesh_.faceOwner()[facei];
328 if (cellRegion[own] != -1)
330 own = cellRegionMaster[cellRegion[own]];
333 label
patchID, zoneID, zoneFlip;
335 getFaceInfo(facei,
patchID, zoneID, zoneFlip);
339 if (mesh_.isInternalFace(facei))
341 nei = mesh_.faceNeighbour()[facei];
343 if (cellRegion[nei] != -1)
345 nei = cellRegionMaster[cellRegion[nei]];
350 DynamicList<label> faceVerts(edgeLoop.size());
354 label pointi = fp.meshPoints()[edgeLoop[i]];
356 if (pointsToRemove.found(pointi))
363 faceVerts.append(pointi);
368 mergedFace.transfer(faceVerts);
405 if (patchFacei != masterIndex)
409 meshMod.setAction(polyRemoveFace(
faceLabels[patchFacei], facei));
416 void Foam::removeFaces::getFaceInfo
427 if (!mesh_.isInternalFace(facei))
429 patchID = mesh_.boundaryMesh().whichPatch(facei);
432 zoneID = mesh_.faceZones().whichZone(facei);
438 const faceZone& fZone = mesh_.faceZones()[zoneID];
440 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
452 const face&
f = mesh_.faces()[facei];
462 if (!pointsToRemove.found(vertI))
464 newFace[newFp++] = vertI;
470 return face(newFace);
475 void Foam::removeFaces::modFace
478 const label masterFaceID,
481 const bool flipFaceFlux,
482 const label newPatchID,
483 const bool removeFromZone,
487 polyTopoChange& meshMod
490 if ((nei == -1) || (own < nei))
562 Foam::removeFaces::removeFaces
589 const labelList& faceOwner = mesh_.faceOwner();
590 const labelList& faceNeighbour = mesh_.faceNeighbour();
592 cellRegion.
setSize(mesh_.nCells());
595 regionMaster.
setSize(mesh_.nCells());
602 label facei = facesToRemove[i];
604 if (!mesh_.isInternalFace(facei))
611 label own = faceOwner[facei];
612 label nei = faceNeighbour[facei];
614 label region0 = cellRegion[own];
615 label region1 = cellRegion[nei];
622 cellRegion[own] = nRegions;
623 cellRegion[nei] = nRegions;
626 regionMaster[nRegions] = own;
632 cellRegion[own] = region1;
634 regionMaster[region1] =
min(own, regionMaster[region1]);
642 cellRegion[nei] = region0;
646 else if (region0 != region1)
649 label freedRegion = -1;
650 label keptRegion = -1;
652 if (region0 < region1)
662 keptRegion = region0;
663 freedRegion = region1;
665 else if (region1 < region0)
675 keptRegion = region1;
676 freedRegion = region0;
679 label master0 = regionMaster[region0];
680 label master1 = regionMaster[region1];
682 regionMaster[freedRegion] = -1;
683 regionMaster[keptRegion] =
min(master0, master1);
688 regionMaster.
setSize(nRegions);
699 label r = cellRegion[celli];
705 if (celli < regionMaster[r])
708 <<
"Not lowest numbered : cell:" << celli
710 <<
" regionmaster:" << regionMaster[r]
718 if (nCells[region] == 1)
721 <<
"Region " << region
722 <<
" has only " << nCells[region] <<
" cells in it" 730 label nUsedRegions = 0;
734 if (regionMaster[i] != -1)
743 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
745 label own = faceOwner[facei];
746 label nei = faceNeighbour[facei];
748 if (cellRegion[own] != -1 && cellRegion[own] == cellRegion[nei])
752 allFacesToRemove.
append(facei);
773 Pout<<
"Writing faces to remove to faceSet " << facesToRemove.
name()
775 facesToRemove.
write();
779 boolList removedFace(mesh_.nFaces(),
false);
785 if (!mesh_.isInternalFace(facei))
788 <<
"Face to remove is not internal face:" << facei
792 removedFace[facei] =
true;
805 labelList faceRegion(mesh_.nFaces(), -1);
811 DynamicList<label> fe;
812 DynamicList<label> ef;
816 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
820 labelList nFacesPerEdge(mesh_.nEdges(), -1);
827 const labelList& fEdges = mesh_.faceEdges(facei, fe);
831 label edgeI = fEdges[i];
833 if (nFacesPerEdge[edgeI] == -1)
835 nFacesPerEdge[edgeI] = mesh_.edgeFaces(edgeI, ef).size()-1;
839 nFacesPerEdge[edgeI]--;
849 forAll(mesh_.edges(), edgeI)
851 if (nFacesPerEdge[edgeI] == -1)
856 const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
858 if (eFaces.size() > 2)
860 nFacesPerEdge[edgeI] = eFaces.
size();
862 else if (eFaces.size() == 2)
868 const edge&
e = mesh_.edges()[edgeI];
871 <<
"Problem : edge has too few face neighbours:" 875 <<
" coords:" << mesh_.points()[
e[0]]
876 << mesh_.points()[
e[1]]
886 OFstream str(mesh_.time().path()/
"edgesWithTwoFaces.obj");
887 Pout<<
"Dumping edgesWithTwoFaces to " << str.
name() <<
endl;
890 forAll(nFacesPerEdge, edgeI)
892 if (nFacesPerEdge[edgeI] == 2)
895 const edge&
e = mesh_.edges()[edgeI];
901 str<<
"l " << vertI-1 <<
' ' << vertI <<
nl;
912 forAll(nFacesPerEdge, edgeI)
914 if (nFacesPerEdge[edgeI] == 2)
920 const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
924 label facei = eFaces[i];
926 if (!removedFace[facei])
940 if (!mesh_.isInternalFace(f0) && !mesh_.isInternalFace(f1))
948 if (patch0 != patch1)
952 <<
"not merging faces " << f0 <<
" and " 953 << f1 <<
" across patch boundary edge " << edgeI
957 nFacesPerEdge[edgeI] = 3;
959 else if (minCos_ < 1 && minCos_ > -1)
961 const polyPatch& pp0 =
patches[patch0];
969 & n0[f1 - pp0.start()]
975 <<
"not merging faces " << f0 <<
" and " 976 << f1 <<
" across edge " << edgeI
981 nFacesPerEdge[edgeI] = 3;
985 else if (mesh_.isInternalFace(f0) != mesh_.isInternalFace(f1))
987 const edge&
e = mesh_.edges()[edgeI];
991 <<
"Problem : edge would have one boundary face" 992 <<
" and one internal face using it." <<
endl 993 <<
"Your remove pattern is probably incorrect." <<
endl 995 <<
" nFaces:" << nFacesPerEdge[edgeI]
997 <<
" coords:" << mesh_.points()[
e[0]]
998 << mesh_.points()[
e[1]]
1014 cellRegion[mesh_.faceOwner()[f0]],
1015 cellRegion[mesh_.faceNeighbour()[f0]]
1019 cellRegion[mesh_.faceOwner()[f1]],
1020 cellRegion[mesh_.faceNeighbour()[f1]]
1023 if (ownEdge != neiEdge)
1025 nFacesPerEdge[edgeI] = 3;
1034 forAll(nFacesPerEdge, edgeI)
1036 if (nFacesPerEdge[edgeI] == 1)
1038 const edge&
e = mesh_.edges()[edgeI];
1041 <<
"Problem : edge would get 1 face using it only" 1042 <<
" edge:" << edgeI
1043 <<
" nFaces:" << nFacesPerEdge[edgeI]
1044 <<
" vertices:" <<
e 1045 <<
" coords:" << mesh_.points()[
e[0]]
1046 <<
' ' << mesh_.points()[
e[1]]
1105 forAll(nFacesPerEdge, edgeI)
1107 if (nFacesPerEdge[edgeI] == 0)
1110 edgesToRemove.insert(edgeI);
1112 else if (nFacesPerEdge[edgeI] == 1)
1116 else if (nFacesPerEdge[edgeI] == 2)
1119 edgesToRemove.insert(edgeI);
1125 OFstream str(mesh_.time().path()/
"edgesToRemove.obj");
1126 Pout<<
"Dumping edgesToRemove to " << str.
name() <<
endl;
1129 for (
const label edgei : edgesToRemove)
1132 const edge&
e = mesh_.edges()[edgei];
1138 str<<
"l " << vertI-1 <<
' ' << vertI <<
nl;
1146 label startFacei = 0;
1151 for (; startFacei < mesh_.nFaces(); startFacei++)
1153 if (faceRegion[startFacei] == -1 && !removedFace[startFacei])
1159 if (startFacei == mesh_.nFaces())
1167 label nRegion = changeFaceRegion
1174 mesh_.faceEdges(startFacei, fe),
1182 else if (nRegion == 1)
1185 faceRegion[startFacei] = -2;
1210 label facei = mesh_.nInternalFaces();
1211 facei < mesh_.nFaces();
1216 label nbrRegion = nbrFaceRegion[facei];
1217 label myRegion = faceRegion[facei];
1219 if (myRegion <= -1 || nbrRegion <= -1)
1221 if (nbrRegion != myRegion)
1224 <<
"Inconsistent face region across coupled patches." 1226 <<
"This side has for facei:" << facei
1227 <<
" region:" << myRegion <<
endl 1228 <<
"The other side has region:" << nbrRegion
1230 <<
"(region -1 means face is to be deleted)" 1234 else if (toNbrRegion[myRegion] == -1)
1237 toNbrRegion[myRegion] = nbrRegion;
1242 if (toNbrRegion[myRegion] != nbrRegion)
1245 <<
"Inconsistent face region across coupled patches." 1247 <<
"This side has for facei:" << facei
1248 <<
" region:" << myRegion
1249 <<
" with coupled neighbouring regions:" 1250 << toNbrRegion[myRegion] <<
" and " 1280 labelList nEdgesPerPoint(mesh_.nPoints());
1284 forAll(pointEdges, pointi)
1286 nEdgesPerPoint[pointi] = pointEdges[pointi].size();
1289 for (
const label edgei : edgesToRemove)
1292 const edge&
e = mesh_.edges()[edgei];
1296 nEdgesPerPoint[
e[i]]--;
1301 forAll(nEdgesPerPoint, pointi)
1303 if (nEdgesPerPoint[pointi] == 1)
1306 <<
"Problem : point would get 1 edge using it only." 1307 <<
" pointi:" << pointi
1308 <<
" coord:" << mesh_.points()[pointi]
1323 forAll(nEdgesPerPoint, pointi)
1325 if (nEdgesPerPoint[pointi] == 0)
1327 pointsToRemove.insert(pointi);
1329 else if (nEdgesPerPoint[pointi] == 1)
1333 else if (nEdgesPerPoint[pointi] == 2)
1336 pointsToRemove.insert(pointi);
1344 OFstream str(mesh_.time().path()/
"pointsToRemove.obj");
1345 Pout<<
"Dumping pointsToRemove to " << str.
name() <<
endl;
1347 for (
const label pointi : pointsToRemove)
1393 if (affectedFace[facei])
1395 affectedFace[facei] =
false;
1397 meshMod.
setAction(polyRemoveFace(facei, -1));
1403 for (
const label pointi : pointsToRemove)
1405 meshMod.
setAction(polyRemovePoint(pointi, -1));
1410 forAll(cellRegion, celli)
1412 label region = cellRegion[celli];
1414 if (region != -1 && (celli != cellRegionMaster[region]))
1416 meshMod.
setAction(polyRemoveCell(celli, cellRegionMaster[region]));
1429 forAll(regionToFaces, regionI)
1431 const labelList& rFaces = regionToFaces[regionI];
1433 if (rFaces.size() <= 1)
1436 <<
"Region:" << regionI
1437 <<
" contains only faces " << rFaces
1453 affectedFace[rFaces[i]] =
false;
1464 forAll(affectedFace, facei)
1466 if (affectedFace[facei])
1468 affectedFace[facei] =
false;
1470 face
f(filterFace(pointsToRemove, facei));
1472 label own = mesh_.faceOwner()[facei];
1474 if (cellRegion[own] != -1)
1476 own = cellRegionMaster[cellRegion[own]];
1479 label
patchID, zoneID, zoneFlip;
1481 getFaceInfo(facei,
patchID, zoneID, zoneFlip);
1485 if (mesh_.isInternalFace(facei))
1487 nei = mesh_.faceNeighbour()[facei];
1489 if (cellRegion[nei] != -1)
1491 nei = cellRegionMaster[cellRegion[nei]];
void size(const label n)
Older name for setAddressableSize.
A face is a list of labels corresponding to mesh vertices.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void transfer(List< T > &list)
Transfer the contents of the argument List into this list and annul the argument list.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
void append(const T &val)
Append an element at the end of the list.
const word & name() const noexcept
Return the object name.
constexpr char nl
The newline '\n' character (0x0a)
Ostream & endl(Ostream &os)
Add newline and flush stream.
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
A PrimitivePatch with an IndirectList for the faces, const reference for the point field...
virtual const fileName & name() const override
Get the name of the output serial stream. (eg, the name of the Fstream file name) ...
labelList faceLabels(nFaceLabels)
List< labelList > labelListList
List of labelList.
#define forAll(list, i)
Loop across all elements in list.
label fcIndex(const label i) const noexcept
The forward circular index. The next index in the list which returns to the first at the end of the l...
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
List< face > faceList
List of faces.
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
void setSize(const label n)
Alias for resize()
labelListList invertOneToMany(const label len, const labelUList &map)
Invert one-to-many map. Unmapped elements will be size 0.
virtual bool write(const bool writeOnProc=true) const
Write using setting from DB.
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
label whichPatch(const label meshFacei) const
Return patch index for a given mesh face index. Uses binary search.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
label find(const T &val) const
Find index of the first occurrence of the value.
int debug
Static debugging option.
defineTypeNameAndDebug(combustionModel, 0)
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells].reset(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
#define WarningInFunction
Report a warning using Foam::Warning.
Direct mesh changes based on v1.3 polyTopoChange syntax.
const polyBoundaryMesh & patches
void setRefinement(const labelList &piercedFaces, const labelList &cellRegion, const labelList &cellRegionMaster, polyTopoChange &) const
Play commands into polyTopoChange to remove faces.
Field< vector > vectorField
Specialisation of Field<T> for vector.
Mesh consisting of general polyhedral cells.
List< label > labelList
A List of labels.
label rcIndex(const label i) const noexcept
The reverse circular index. The previous index in the list which returns to the last at the beginning...
List< bool > boolList
A List of bools.
const labelListList & cellCells() const
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
label compatibleRemoves(const labelList &inPiercedFaces, labelList &cellRegion, labelList &cellRegionMaster, labelList &outPiercedFaces) const
Find faces including those with cells which have the same mastercell.
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.
static constexpr const zero Zero
Global zero (0)