52 if (
findIndex(elems2, elems1[elemI]) != -1)
61 bool Foam::meshCutter::isIn
76 cuts[cuts.fcIndex(index)] == twoCuts[1]
77 || cuts[cuts.rcIndex(index)] == twoCuts[1]
94 if (cuts.cellLoops()[celli].size())
103 Foam::label Foam::meshCutter::findInternalFacePoint
118 if (mesh().isInternalFace(facei))
135 void Foam::meshCutter::faceCells
146 const face&
f = mesh().faces()[facei];
148 own = mesh().faceOwner()[facei];
150 if (cellLoops[own].size() && uses(
f, anchorPts[own]))
152 own = addedCells_[own];
157 if (mesh().isInternalFace(facei))
159 nei = mesh().faceNeighbour()[facei];
161 if (cellLoops[nei].size() && uses(
f, anchorPts[nei]))
163 nei = addedCells_[nei];
169 void Foam::meshCutter::getFaceInfo
179 if (!mesh().isInternalFace(facei))
181 patchID = mesh().boundaryMesh().whichPatch(facei);
184 zoneID = mesh().faceZones().whichZone(facei);
190 const faceZone& fZone = mesh().faceZones()[zoneID];
192 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
197 void Foam::meshCutter::addFace
199 polyTopoChange& meshMod,
206 label patchID, zoneID, zoneFlip;
208 getFaceInfo(facei, patchID, zoneID, zoneFlip);
210 if ((nei == -1) || (own < nei))
215 Pout<<
"Adding face " << newFace
216 <<
" with new owner:" << own
217 <<
" with new neighbour:" << nei
218 <<
" patchID:" << patchID
219 <<
" zoneID:" << zoneID
220 <<
" zoneFlip:" << zoneFlip
246 Pout<<
"Adding (reversed) face " << newFace.reverseFace()
247 <<
" with new owner:" << nei
248 <<
" with new neighbour:" << own
249 <<
" patchID:" << patchID
250 <<
" zoneID:" << zoneID
251 <<
" zoneFlip:" << zoneFlip
259 newFace.reverseFace(),
275 void Foam::meshCutter::modFace
277 polyTopoChange& meshMod,
284 label patchID, zoneID, zoneFlip;
286 getFaceInfo(facei, patchID, zoneID, zoneFlip);
290 (own != mesh().faceOwner()[facei])
292 mesh().isInternalFace(facei)
293 && (nei != mesh().faceNeighbour()[facei])
295 || (newFace != mesh().faces()[facei])
300 Pout<<
"Modifying face " << facei
301 <<
" old vertices:" << mesh().faces()[facei]
302 <<
" new vertices:" << newFace
303 <<
" new owner:" << own
304 <<
" new neighbour:" << nei
305 <<
" new zoneID:" << zoneID
306 <<
" new zoneFlip:" << zoneFlip
310 if ((nei == -1) || (own < nei))
334 newFace.reverseFace(),
350 void Foam::meshCutter::copyFace
364 newFace[newFp++] =
f[fp];
366 fp = (fp + 1) %
f.
size();
368 newFace[newFp] =
f[fp];
372 void Foam::meshCutter::splitFace
388 <<
"Cannot find vertex (new numbering) " << v0
398 <<
"Cannot find vertex (new numbering) " << v1
404 f0.setSize((endFp + 1 +
f.
size() - startFp) %
f.
size());
405 f1.setSize(
f.
size() - f0.size() + 2);
407 copyFace(
f, startFp, endFp, f0);
408 copyFace(
f, endFp, startFp, f1);
412 Foam::face Foam::meshCutter::addEdgeCutsToFace(
const label facei)
const
414 const face&
f = mesh().faces()[facei];
416 face newFace(2 *
f.
size());
423 newFace[newFp++] =
f[fp];
428 HashTable<label, edge, Hash<edge>>::const_iterator fnd =
429 addedPoints_.find(edge(
f[fp],
f[fp1]));
431 if (fnd != addedPoints_.end())
434 newFace[newFp++] = fnd();
438 newFace.setSize(newFp);
450 face newFace(2*loop.size());
456 label cut = loop[fp];
460 label edgeI = getEdge(cut);
462 const edge&
e = mesh().edges()[edgeI];
464 label vertI = addedPoints_[
e];
466 newFace[newFacei++] = vertI;
471 label vertI = getVertex(cut);
473 newFace[newFacei++] = vertI;
475 label nextCut = loop[loop.fcIndex(fp)];
477 if (!isEdge(nextCut))
480 label nextVertI = getVertex(nextCut);
487 HashTable<label, edge, Hash<edge>>::const_iterator fnd =
488 addedPoints_.find(mesh().edges()[edgeI]);
490 if (fnd != addedPoints_.end())
492 newFace[newFacei++] = fnd();
498 newFace.setSize(newFacei);
532 addedCells_.resize(cuts.
nLoops());
535 addedFaces_.resize(cuts.
nLoops());
537 addedPoints_.clear();
538 addedPoints_.resize(cuts.
nLoops());
551 boolList edgeOnCutCell(mesh().nEdges(),
false);
556 const labelList& cEdges = mesh().cellEdges(celli);
559 edgeOnCutCell[cEdges[i]] =
true;
567 if (cuts.
edgeIsCut()[edgeI] && !edgeOnCutCell[edgeI])
569 const edge&
e = mesh().edges()[edgeI];
572 <<
"Problem: cut edge but none of the cells using"
574 <<
"edge:" << edgeI <<
" verts:" <<
e
575 <<
" at:" <<
e.line(mesh().
points())
590 const edge&
e = mesh().edges()[edgeI];
593 label masterPointi =
e.start();
595 const point& v0 = mesh().points()[
e.start()];
596 const point& v1 = mesh().points()[
e.end()];
600 point newPt = weight*v1 + (1.0-weight)*v0;
615 addedPoints_.insert(
e, addedPointi);
619 Pout<<
"Added point " << addedPointi
621 << masterPointi <<
" of edge " << edgeI
622 <<
" vertices " <<
e <<
endl;
633 if (cellLoops[celli].size())
645 mesh().cellZones().whichZone(celli)
649 addedCells_.insert(celli, addedCelli);
653 Pout<<
"Added cell " << addedCells_[celli] <<
" to cell "
666 const labelList& loop = cellLoops[celli];
673 face newFace(loopToFace(celli, loop));
676 label masterPointi = findInternalFacePoint(anchorPts[celli]);
696 addedFaces_.insert(celli, addedFacei);
714 Pout<<
"Added splitting face " << newFace <<
" index:"
716 <<
" to owner " << celli
717 <<
" neighbour " << addedCells_[celli]
719 writeCuts(
Pout, loop, weights);
733 boolList faceUptodate(mesh().nFaces(),
false);
739 label facei = iter.key();
742 face newFace(addEdgeCutsToFace(facei));
745 const edge& splitEdge = iter();
747 label cut0 = splitEdge[0];
752 label edgeI = getEdge(cut0);
753 v0 = addedPoints_[mesh().edges()[edgeI]];
757 v0 = getVertex(cut0);
760 label cut1 = splitEdge[1];
764 label edgeI = getEdge(cut1);
765 v1 = addedPoints_[mesh().edges()[edgeI]];
769 v1 = getVertex(cut1);
774 splitFace(newFace, v0, v1, f0, f1);
776 label own = mesh().faceOwner()[facei];
780 if (mesh().isInternalFace(facei))
782 nei = mesh().faceNeighbour()[facei];
787 Pout<<
"Split face " << mesh().faces()[facei]
788 <<
" own:" << own <<
" nei:" << nei
790 <<
" and f1:" << f1 <<
endl;
801 const face&
f = mesh().faces()[facei];
806 if (cellLoops[own].empty())
811 else if (isIn(splitEdge, cellLoops[own]))
815 if (uses(f0, anchorPts[own]))
817 f0Owner = addedCells_[own];
823 f1Owner = addedCells_[own];
830 if (uses(
f, anchorPts[own]))
832 label newCelli = addedCells_[own];
844 label f0Neighbour = -1;
845 label f1Neighbour = -1;
849 if (cellLoops[nei].empty())
854 else if (isIn(splitEdge, cellLoops[nei]))
858 if (uses(f0, anchorPts[nei]))
860 f0Neighbour = addedCells_[nei];
866 f1Neighbour = addedCells_[nei];
873 if (uses(
f, anchorPts[nei]))
875 label newCelli = addedCells_[nei];
876 f0Neighbour = newCelli;
877 f1Neighbour = newCelli;
888 addFace(meshMod, facei, f0, f0Owner, f0Neighbour);
890 modFace(meshMod, facei, f1, f1Owner, f1Neighbour);
892 faceUptodate[facei] =
true;
905 if (edgeIsCut[edgeI])
907 const labelList& eFaces = mesh().edgeFaces()[edgeI];
911 label facei = eFaces[i];
913 if (!faceUptodate[facei])
916 face newFace(addEdgeCutsToFace(facei));
920 Pout<<
"Added edge cuts to face " << facei
921 <<
" f:" << mesh().faces()[facei]
922 <<
" newFace:" << newFace <<
endl;
927 faceCells(cuts, facei, own, nei);
929 modFace(meshMod, facei, newFace, own, nei);
931 faceUptodate[facei] =
true;
944 if (cellLoops[celli].size())
946 const labelList& cllFaces = mesh().cells()[celli];
948 forAll(cllFaces, cllFacei)
950 label facei = cllFaces[cllFacei];
952 if (!faceUptodate[facei])
955 const face&
f = mesh().faces()[facei];
957 if (debug && (
f != addEdgeCutsToFace(facei)))
960 <<
"Problem: edges added to face which does "
961 <<
" not use a marked cut" <<
endl
962 <<
"facei:" << facei <<
endl
963 <<
"face:" <<
f <<
endl
964 <<
"newFace:" << addEdgeCutsToFace(facei)
970 faceCells(cuts, facei, own, nei);
981 faceUptodate[facei] =
true;
989 Pout<<
"meshCutter:" <<
nl
990 <<
" cells split:" << addedCells_.size() <<
nl
991 <<
" faces added:" << addedFaces_.size() <<
nl
992 <<
" points added on edges:" << addedPoints_.size() <<
nl
1005 Map<label> newAddedCells(addedCells_.size());
1009 label celli = iter.key();
1012 label addedCelli = iter();
1016 if (newCelli >= 0 && newAddedCelli >= 0)
1021 && (newCelli != celli || newAddedCelli != addedCelli)
1024 Pout<<
"meshCutter::topoChange :"
1025 <<
" updating addedCell for cell " << celli
1026 <<
" from " << addedCelli
1027 <<
" to " << newAddedCelli <<
endl;
1029 newAddedCells.
insert(newCelli, newAddedCelli);
1034 addedCells_.transfer(newAddedCells);
1038 Map<label> newAddedFaces(addedFaces_.size());
1042 label celli = iter.key();
1045 label addedFacei = iter();
1049 if ((newCelli >= 0) && (newAddedFacei >= 0))
1054 && (newCelli != celli || newAddedFacei != addedFacei)
1057 Pout<<
"meshCutter::topoChange :"
1058 <<
" updating addedFace for cell " << celli
1059 <<
" from " << addedFacei
1060 <<
" to " << newAddedFacei
1063 newAddedFaces.
insert(newCelli, newAddedFacei);
1068 addedFaces_.transfer(newAddedFaces);
1077 addedPoints_.begin();
1078 iter != addedPoints_.end();
1082 const edge&
e = iter.key();
1088 label addedPointi = iter();
1092 if ((newStart >= 0) && (newEnd >= 0) && (newAddedPointi >= 0))
1094 edge newE =
edge(newStart, newEnd);
1099 && (
e != newE || newAddedPointi != addedPointi)
1102 Pout<<
"meshCutter::topoChange :"
1103 <<
" updating addedPoints for edge " <<
e
1104 <<
" from " << addedPointi
1105 <<
" to " << newAddedPointi
1109 newAddedPoints.
insert(newE, newAddedPointi);
1114 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)
Description of cuts across cells.
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.
meshCutter(const polyMesh &mesh)
Construct from mesh.
void topoChange(const polyTopoChangeMap &)
Force recalculation of locally stored data on topological change.
void setRefinement(const cellCuts &cuts, polyTopoChange &meshMod)
Do actual cutting with cut description. Inserts mesh changes.
Class containing data for cell addition.
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 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.
#define WarningInFunction
Report a warning using Foam::Warning.
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)
List< labelList > labelListList
A List of labelList.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
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]