69 bool Foam::meshCutAndRemove::isIn
84 cuts[cuts.fcIndex(index)] == twoCuts[1]
85 || cuts[cuts.rcIndex(index)] == twoCuts[1]
102 if (cuts.cellLoops()[celli].size())
111 Foam::label Foam::meshCutAndRemove::findInternalFacePoint
126 if (mesh().isInternalFace(facei))
143 Foam::label Foam::meshCutAndRemove::findPatchFacePoint
146 const label exposedPatchi
150 const polyBoundaryMesh&
patches = mesh().boundaryMesh();
173 void Foam::meshCutAndRemove::faceCells
176 const label exposedPatchi,
186 const face&
f = mesh().faces()[facei];
188 own = mesh().faceOwner()[facei];
190 if (cellLoops[own].size() && firstCommon(
f, anchorPts[own]) == -1)
198 if (mesh().isInternalFace(facei))
200 nei = mesh().faceNeighbour()[facei];
202 if (cellLoops[nei].size() && firstCommon(
f, anchorPts[nei]) == -1)
208 patchID = mesh().boundaryMesh().whichPatch(facei);
210 if (patchID == -1 && (own == -1 || nei == -1))
213 patchID = exposedPatchi;
218 void Foam::meshCutAndRemove::getZoneInfo
225 zoneID = mesh().faceZones().whichZone(facei);
231 const faceZone& fZone = mesh().faceZones()[zoneID];
233 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
238 void Foam::meshCutAndRemove::addFace
240 polyTopoChange& meshMod,
242 const label masterPointi,
252 getZoneInfo(facei, zoneID, zoneFlip);
254 if ((nei == -1) || (own != -1 && own < nei))
259 Pout<<
"Adding face " << newFace
260 <<
" with new owner:" << own
261 <<
" with new neighbour:" << nei
262 <<
" patchID:" << patchID
263 <<
" anchor:" << masterPointi
264 <<
" zoneID:" << zoneID
265 <<
" zoneFlip:" << zoneFlip
291 Pout<<
"Adding (reversed) face " << newFace.reverseFace()
292 <<
" with new owner:" << nei
293 <<
" with new neighbour:" << own
294 <<
" patchID:" << patchID
295 <<
" anchor:" << masterPointi
296 <<
" zoneID:" << zoneID
297 <<
" zoneFlip:" << zoneFlip
305 newFace.reverseFace(),
322 void Foam::meshCutAndRemove::modFace
324 polyTopoChange& meshMod,
335 getZoneInfo(facei, zoneID, zoneFlip);
339 (own != mesh().faceOwner()[facei])
341 mesh().isInternalFace(facei)
342 && (nei != mesh().faceNeighbour()[facei])
344 || (newFace != mesh().faces()[facei])
349 Pout<<
"Modifying face " << facei
350 <<
" old vertices:" << mesh().faces()[facei]
351 <<
" new vertices:" << newFace
352 <<
" new owner:" << own
353 <<
" new neighbour:" << nei
354 <<
" new patch:" << patchID
355 <<
" new zoneID:" << zoneID
356 <<
" new zoneFlip:" << zoneFlip
360 if ((nei == -1) || (own != -1 && own < nei))
384 newFace.reverseFace(),
400 void Foam::meshCutAndRemove::copyFace
414 newFace[newFp++] =
f[fp];
416 fp = (fp + 1) %
f.
size();
418 newFace[newFp] =
f[fp];
426 void Foam::meshCutAndRemove::splitFace
442 <<
"Cannot find vertex (new numbering) " << v0
452 <<
"Cannot find vertex (new numbering) " << v1
458 f0.setSize((endFp + 1 +
f.
size() - startFp) %
f.
size());
459 f1.setSize(
f.
size() - f0.size() + 2);
461 copyFace(
f, startFp, endFp, f0);
462 copyFace(
f, endFp, startFp, f1);
466 Foam::face Foam::meshCutAndRemove::addEdgeCutsToFace(
const label facei)
const
468 const face&
f = mesh().faces()[facei];
470 face newFace(2 *
f.
size());
477 newFace[newFp++] =
f[fp];
482 HashTable<label, edge, Hash<edge>>::const_iterator fnd =
483 addedPoints_.find(edge(
f[fp],
f[fp1]));
485 if (fnd != addedPoints_.end())
488 newFace[newFp++] = fnd();
492 newFace.setSize(newFp);
507 face newFace(2*loop.size());
513 label cut = loop[fp];
517 label edgeI = getEdge(cut);
519 const edge&
e = mesh().edges()[edgeI];
521 label vertI = addedPoints_[
e];
523 newFace[newFacei++] = vertI;
528 label vertI = getVertex(cut);
530 newFace[newFacei++] = vertI;
532 label nextCut = loop[loop.fcIndex(fp)];
534 if (!isEdge(nextCut))
537 label nextVertI = getVertex(nextCut);
544 HashTable<label, edge, Hash<edge>>::const_iterator fnd =
545 addedPoints_.find(mesh().edges()[edgeI]);
547 if (fnd != addedPoints_.end())
549 newFace[newFacei++] = fnd();
555 newFace.setSize(newFacei);
576 const label exposedPatchi,
584 addedFaces_.resize(cuts.
nLoops());
586 addedPoints_.clear();
587 addedPoints_.resize(cuts.
nLoops());
598 if (exposedPatchi < 0 || exposedPatchi >=
patches.
size())
601 <<
"Illegal exposed patch " << exposedPatchi
614 const edge&
e = mesh().edges()[edgeI];
617 if (debug && findCutCell(cuts, mesh().edgeCells()[edgeI]) == -1)
620 <<
"Problem: cut edge but none of the cells using it is\n"
621 <<
"edge:" << edgeI <<
" verts:" <<
e
626 label masterPointi =
e.start();
628 const point& v0 = mesh().points()[
e.start()];
629 const point& v1 = mesh().points()[
e.end()];
633 point newPt = weight*v1 + (1.0-weight)*v0;
648 addedPoints_.insert(
e, addedPointi);
652 Pout<<
"Added point " << addedPointi
654 << masterPointi <<
" of edge " << edgeI
655 <<
" vertices " <<
e <<
endl;
669 const labelList& loop = cellLoops[celli];
676 label cut = loop[fp];
680 usedPoint[getVertex(cut)] =
true;
684 const labelList& anchors = anchorPts[celli];
688 usedPoint[anchors[i]] =
true;
694 const labelList& cPoints = mesh().cellPoints()[celli];
698 usedPoint[cPoints[i]] =
true;
709 const edge& fCut = iter();
717 label pointi = getVertex(cut);
719 if (!usedPoint[pointi])
722 <<
"Problem: faceSplitCut not used by any loop"
723 <<
" or cell anchor point"
724 <<
"face:" << iter.key() <<
" point:" << pointi
725 <<
" coord:" << mesh().points()[pointi]
736 if (!usedPoint[pointi])
739 <<
"Problem: point is marked as cut but"
740 <<
" not used by any loop"
741 <<
" or cell anchor point"
742 <<
"point:" << pointi
743 <<
" coord:" << mesh().points()[pointi]
753 if (!usedPoint[pointi])
759 Pout<<
"Removing unused point " << pointi <<
endl;
772 const labelList& loop = cellLoops[celli];
776 if (cutPatch[celli] < 0 || cutPatch[celli] >=
patches.
size())
779 <<
"Illegal patch " << cutPatch[celli]
780 <<
" provided for cut cell " << celli
788 face newFace(loopToFace(celli, loop));
793 label masterPointi = findPatchFacePoint(newFace, exposedPatchi);
813 addedFaces_.insert(celli, addedFacei);
817 Pout<<
"Added splitting face " << newFace <<
" index:"
818 << addedFacei <<
" from masterPoint:" << masterPointi
819 <<
" to owner " << celli <<
" with anchors:"
836 writeCuts(
Pout, loop, weights);
851 boolList faceUptodate(mesh().nFaces(),
false);
857 label facei = iter.key();
860 face newFace(addEdgeCutsToFace(facei));
863 const edge& splitEdge = iter();
865 label cut0 = splitEdge[0];
870 label edgeI = getEdge(cut0);
871 v0 = addedPoints_[mesh().edges()[edgeI]];
875 v0 = getVertex(cut0);
878 label cut1 = splitEdge[1];
882 label edgeI = getEdge(cut1);
883 v1 = addedPoints_[mesh().edges()[edgeI]];
887 v1 = getVertex(cut1);
892 splitFace(newFace, v0, v1, f0, f1);
894 label own = mesh().faceOwner()[facei];
898 if (mesh().isInternalFace(facei))
900 nei = mesh().faceNeighbour()[facei];
905 Pout<<
"Split face " << mesh().faces()[facei]
906 <<
" own:" << own <<
" nei:" << nei
908 <<
" and f1:" << f1 <<
endl;
922 const face&
f = mesh().faces()[facei];
927 if (cellLoops[own].empty())
933 else if (isIn(splitEdge, cellLoops[own]))
937 if (firstCommon(f0, anchorPts[own]) != -1)
954 if (firstCommon(
f, anchorPts[own]) != -1)
974 if (cellLoops[nei].empty())
979 else if (isIn(splitEdge, cellLoops[nei]))
985 if (firstCommon(f0, anchorPts[nei]) != -1)
1001 if (firstCommon(
f, anchorPts[nei]) != -1)
1018 Pout<<
"f0 own:" << f0Own <<
" nei:" << f0Nei
1019 <<
" f1 own:" << f1Own <<
" nei:" << f1Nei
1030 patchID = exposedPatchi;
1037 bool modifiedFacei =
false;
1044 modFace(meshMod, facei, f0, f0Own, f0Nei, patchID);
1045 modifiedFacei =
true;
1053 modFace(meshMod, facei, f0, f0Own, f0Nei, patchID);
1054 modifiedFacei =
true;
1059 modFace(meshMod, facei, f0, f0Own, f0Nei, -1);
1060 modifiedFacei =
true;
1078 modFace(meshMod, facei, f1, f1Own, f1Nei, patchID);
1079 modifiedFacei =
true;
1083 label masterPointi = findPatchFacePoint(f1, patchID);
1105 modFace(meshMod, facei, f1, f1Own, f1Nei, patchID);
1106 modifiedFacei =
true;
1110 label masterPointi = findPatchFacePoint(f1, patchID);
1129 modFace(meshMod, facei, f1, f1Own, f1Nei, -1);
1130 modifiedFacei =
true;
1134 label masterPointi = findPatchFacePoint(f1, -1);
1136 addFace(meshMod, facei, masterPointi, f1, f1Own, f1Nei, -1);
1141 if (f0Own == -1 && f0Nei == -1 && !modifiedFacei)
1147 Pout<<
"Removed face " << facei <<
endl;
1151 faceUptodate[facei] =
true;
1164 if (edgeIsCut[edgeI])
1166 const labelList& eFaces = mesh().edgeFaces()[edgeI];
1170 label facei = eFaces[i];
1172 if (!faceUptodate[facei])
1179 label own, nei, patchID;
1180 faceCells(cuts, exposedPatchi, facei, own, nei, patchID);
1183 if (own == -1 && nei == -1)
1189 Pout<<
"Removed face " << facei <<
endl;
1195 face newFace(addEdgeCutsToFace(facei));
1199 Pout<<
"Added edge cuts to face " << facei
1200 <<
" f:" << mesh().faces()[facei]
1201 <<
" newFace:" << newFace <<
endl;
1215 faceUptodate[facei] =
true;
1228 const faceList& faces = mesh().faces();
1232 if (!faceUptodate[facei])
1235 label own, nei, patchID;
1236 faceCells(cuts, exposedPatchi, facei, own, nei, patchID);
1238 if (own == -1 && nei == -1)
1244 Pout<<
"Removed face " << facei <<
endl;
1249 modFace(meshMod, facei, faces[facei], own, nei, patchID);
1252 faceUptodate[facei] =
true;
1258 Pout<<
"meshCutAndRemove:" <<
nl
1259 <<
" cells split:" << cuts.
nLoops() <<
nl
1260 <<
" faces added:" << addedFaces_.size() <<
nl
1261 <<
" points added on edges:" << addedPoints_.size() <<
nl
1271 Map<label> newAddedFaces(addedFaces_.size());
1275 label celli = iter.key();
1278 label addedFacei = iter();
1282 if ((newCelli >= 0) && (newAddedFacei >= 0))
1287 && (newCelli != celli || newAddedFacei != addedFacei)
1290 Pout<<
"meshCutAndRemove::topoChange :"
1291 <<
" updating addedFace for cell " << celli
1292 <<
" from " << addedFacei
1293 <<
" to " << newAddedFacei
1296 newAddedFaces.
insert(newCelli, newAddedFacei);
1301 addedFaces_.transfer(newAddedFaces);
1310 addedPoints_.begin();
1311 iter != addedPoints_.end();
1315 const edge&
e = iter.key();
1321 label addedPointi = iter();
1325 if ((newStart >= 0) && (newEnd >= 0) && (newAddedPointi >= 0))
1327 edge newE =
edge(newStart, newEnd);
1332 && (
e != newE || newAddedPointi != addedPointi)
1335 Pout<<
"meshCutAndRemove::topoChange :"
1336 <<
" updating addedPoints for edge " <<
e
1337 <<
" from " << addedPointi
1338 <<
" to " << newAddedPointi
1342 newAddedPoints.
insert(newE, newAddedPointi);
1347 addedPoints_.transfer(newAddedPoints);
#define forAll(list, i)
Loop across all elements in list.
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
An STL-conforming hash table.
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
Hash function class for primitives. All non-primitives used to hash entries on hash tables likely nee...
void size(const label)
Override size to be inconsistent with allocated storage.
A HashTable to objects of type <T> with a label key.
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
bool empty() const
Return true if the UList is empty (ie, size() is zero)
label size() const
Return the number of elements in the UPtrList.
Description of cuts across cells.
const boolList & pointIsCut() const
Is mesh point cut.
const scalarField & edgeWeight() const
If edge is cut gives weight (ratio between start() and end())
const labelListList & cellLoops() const
For each cut cell the cut along the circumference.
const Map< edge > & faceSplitCut() const
Gives for split face the two cuts that split the face into two.
label nLoops() const
Number of valid cell loops.
const boolList & edgeIsCut() const
Is edge cut.
const labelListList & cellAnchorPoints() const
For each cut cell the points on the 'anchor' side of the cell.
Combines edge or vertex in single label. Used to specify cuts across cell circumference.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
A face is a list of labels corresponding to mesh vertices.
like meshCutter but also removes non-anchor side of cell.
void topoChange(const polyTopoChangeMap &)
Force recalculation of locally stored data on topological change.
void setRefinement(const label exposedPatchi, const cellCuts &cuts, const labelList &cutPatch, polyTopoChange &meshMod)
Do actual cutting with cut description. Inserts mesh changes.
meshCutAndRemove(const polyMesh &mesh)
Construct from mesh.
A face addition data class. A face can be inflated either from a point or from another face and can e...
Class containing data for point addition.
Mesh consisting of general polyhedral cells.
Class containing data for face removal.
Class containing data for point removal.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
const labelList & reversePointMap() const
Reverse point map.
const labelList & reverseCellMap() const
Reverse cell map.
const labelList & reverseFaceMap() const
Reverse face map.
Direct mesh changes based on v1.3 polyTopoChange syntax.
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const fvPatchList & patches
labelListList cellCuts(const cell &c, const cellEdgeAddressing &cAddr, const faceUList &fs, const List< List< labelPair >> &fCuts, const scalarField &pAlphas, const scalar isoAlpha)
Return the cuts for a given cell. This returns a list of lists of cell-edge.
List< label > labelList
A List of labels.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Ostream & endl(Ostream &os)
Add newline and flush stream.
static const labelSphericalTensor labelI(1)
Identity labelTensor.
errorManip< error > abort(error &err)
void reverse(UList< T > &, const label n)
List< labelList > labelListList
A List of labelList.
defineTypeNameAndDebug(combustionModel, 0)
prefixOSstream Pout(cout, "Pout")
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
labelList pointLabels(nPoints, -1)
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=*(cellModeller::lookup("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]=cellShape(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, &mergedCyclicPolyPatch::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]