48 if (
findIndex(elems2, elems1[elemI]) != -1)
57 bool Foam::meshCutter::isIn
72 cuts[cuts.fcIndex(index)] == twoCuts[1]
73 || cuts[cuts.rcIndex(index)] == twoCuts[1]
90 if (cuts.cellLoops()[celli].size())
114 if (mesh().isInternalFace(facei))
131 void Foam::meshCutter::faceCells
142 const face&
f = mesh().faces()[facei];
144 own = mesh().faceOwner()[facei];
146 if (cellLoops[own].size() && uses(
f, anchorPts[own]))
148 own = addedCells_[own];
153 if (mesh().isInternalFace(facei))
155 nei = mesh().faceNeighbour()[facei];
157 if (cellLoops[nei].size() && uses(
f, anchorPts[nei]))
159 nei = addedCells_[nei];
169 if (!mesh().isInternalFace(facei))
171 patchID = mesh().boundaryMesh().whichPatch(facei);
178 void Foam::meshCutter::addFace
180 polyTopoChange& meshMod,
187 const label patchID = getPatchIndex(facei);
189 if ((nei == -1) || (own < nei))
194 Pout<<
"Adding face " << newFace
195 <<
" with new owner:" << own
196 <<
" with new neighbour:" << nei
197 <<
" patchID:" << patchID
216 Pout<<
"Adding (reversed) face " << newFace.reverseFace()
217 <<
" with new owner:" << nei
218 <<
" with new neighbour:" << own
219 <<
" patchID:" << patchID
225 newFace.reverseFace(),
236 void Foam::meshCutter::modifyFace
238 polyTopoChange& meshMod,
245 const label patchID = getPatchIndex(facei);
249 (own != mesh().faceOwner()[facei])
251 mesh().isInternalFace(facei)
252 && (nei != mesh().faceNeighbour()[facei])
254 || (newFace != mesh().faces()[facei])
259 Pout<<
"Modifying face " << facei
260 <<
" old vertices:" << mesh().faces()[facei]
261 <<
" new vertices:" << newFace
262 <<
" new owner:" << own
263 <<
" new neighbour:" << nei
267 if ((nei == -1) || (own < nei))
283 newFace.reverseFace(),
295 void Foam::meshCutter::copyFace
309 newFace[newFp++] =
f[fp];
311 fp = (fp + 1) %
f.
size();
313 newFace[newFp] =
f[fp];
317 void Foam::meshCutter::splitFace
333 <<
"Cannot find vertex (new numbering) " << v0
343 <<
"Cannot find vertex (new numbering) " << v1
349 f0.setSize((endFp + 1 +
f.
size() - startFp) %
f.
size());
350 f1.setSize(
f.
size() - f0.size() + 2);
352 copyFace(
f, startFp, endFp, f0);
353 copyFace(
f, endFp, startFp, f1);
357 Foam::face Foam::meshCutter::addEdgeCutsToFace(
const label facei)
const
359 const face&
f = mesh().faces()[facei];
361 face newFace(2 *
f.
size());
368 newFace[newFp++] =
f[fp];
373 HashTable<label, edge, Hash<edge>>::const_iterator fnd =
374 addedPoints_.find(edge(
f[fp],
f[fp1]));
376 if (fnd != addedPoints_.end())
379 newFace[newFp++] = fnd();
383 newFace.setSize(newFp);
395 face newFace(2*loop.size());
401 label cut = loop[fp];
405 label edgeI = getEdge(cut);
407 const edge&
e = mesh().edges()[edgeI];
409 label vertI = addedPoints_[
e];
411 newFace[newFacei++] = vertI;
416 label vertI = getVertex(cut);
418 newFace[newFacei++] = vertI;
420 label nextCut = loop[loop.fcIndex(fp)];
422 if (!isEdge(nextCut))
425 label nextVertI = getVertex(nextCut);
432 HashTable<label, edge, Hash<edge>>::const_iterator fnd =
433 addedPoints_.find(mesh().edges()[edgeI]);
435 if (fnd != addedPoints_.end())
437 newFace[newFacei++] = fnd();
443 newFace.setSize(newFacei);
477 addedCells_.resize(cuts.
nLoops());
480 addedFaces_.resize(cuts.
nLoops());
482 addedPoints_.clear();
483 addedPoints_.resize(cuts.
nLoops());
496 boolList edgeOnCutCell(mesh().nEdges(),
false);
501 const labelList& cEdges = mesh().cellEdges(celli);
504 edgeOnCutCell[cEdges[i]] =
true;
512 if (cuts.
edgeIsCut()[edgeI] && !edgeOnCutCell[edgeI])
514 const edge&
e = mesh().edges()[edgeI];
517 <<
"Problem: cut edge but none of the cells using"
519 <<
"edge:" << edgeI <<
" verts:" <<
e
520 <<
" at:" <<
e.line(mesh().
points())
535 const edge&
e = mesh().edges()[edgeI];
538 label masterPointi =
e.start();
540 const point& v0 = mesh().points()[
e.start()];
541 const point& v1 = mesh().points()[
e.end()];
545 point newPt = weight*v1 + (1.0-weight)*v0;
555 addedPoints_.insert(
e, addedPointi);
559 Pout<<
"Added point " << addedPointi
561 << masterPointi <<
" of edge " << edgeI
562 <<
" vertices " <<
e <<
endl;
573 if (cellLoops[celli].size())
578 addedCells_.insert(celli, addedCelli);
582 Pout<<
"Added cell " << addedCells_[celli] <<
" to cell "
595 const labelList& loop = cellLoops[celli];
602 const face newFace(loopToFace(celli, loop));
614 addedFaces_.insert(celli, addedFacei);
632 Pout<<
"Added splitting face " << newFace <<
" index:"
634 <<
" to owner " << celli
635 <<
" neighbour " << addedCells_[celli]
637 writeCuts(
Pout, loop, weights);
651 boolList faceUptodate(mesh().nFaces(),
false);
657 label facei = iter.key();
660 face newFace(addEdgeCutsToFace(facei));
663 const edge& splitEdge = iter();
665 label cut0 = splitEdge[0];
670 label edgeI = getEdge(cut0);
671 v0 = addedPoints_[mesh().edges()[edgeI]];
675 v0 = getVertex(cut0);
678 label cut1 = splitEdge[1];
682 label edgeI = getEdge(cut1);
683 v1 = addedPoints_[mesh().edges()[edgeI]];
687 v1 = getVertex(cut1);
692 splitFace(newFace, v0, v1, f0, f1);
694 label own = mesh().faceOwner()[facei];
698 if (mesh().isInternalFace(facei))
700 nei = mesh().faceNeighbour()[facei];
705 Pout<<
"Split face " << mesh().faces()[facei]
706 <<
" own:" << own <<
" nei:" << nei
708 <<
" and f1:" << f1 <<
endl;
719 const face&
f = mesh().faces()[facei];
724 if (cellLoops[own].empty())
729 else if (isIn(splitEdge, cellLoops[own]))
733 if (uses(f0, anchorPts[own]))
735 f0Owner = addedCells_[own];
741 f1Owner = addedCells_[own];
748 if (uses(
f, anchorPts[own]))
750 label newCelli = addedCells_[own];
762 label f0Neighbour = -1;
763 label f1Neighbour = -1;
767 if (cellLoops[nei].empty())
772 else if (isIn(splitEdge, cellLoops[nei]))
776 if (uses(f0, anchorPts[nei]))
778 f0Neighbour = addedCells_[nei];
784 f1Neighbour = addedCells_[nei];
791 if (uses(
f, anchorPts[nei]))
793 label newCelli = addedCells_[nei];
794 f0Neighbour = newCelli;
795 f1Neighbour = newCelli;
806 addFace(meshMod, facei, f0, f0Owner, f0Neighbour);
808 modifyFace(meshMod, facei, f1, f1Owner, f1Neighbour);
810 faceUptodate[facei] =
true;
823 if (edgeIsCut[edgeI])
825 const labelList& eFaces = mesh().edgeFaces()[edgeI];
829 label facei = eFaces[i];
831 if (!faceUptodate[facei])
834 face newFace(addEdgeCutsToFace(facei));
838 Pout<<
"Added edge cuts to face " << facei
839 <<
" f:" << mesh().faces()[facei]
840 <<
" newFace:" << newFace <<
endl;
845 faceCells(cuts, facei, own, nei);
847 modifyFace(meshMod, facei, newFace, own, nei);
849 faceUptodate[facei] =
true;
862 if (cellLoops[celli].size())
864 const labelList& cllFaces = mesh().cells()[celli];
866 forAll(cllFaces, cllFacei)
868 label facei = cllFaces[cllFacei];
870 if (!faceUptodate[facei])
873 const face&
f = mesh().faces()[facei];
875 if (debug && (
f != addEdgeCutsToFace(facei)))
878 <<
"Problem: edges added to face which does "
879 <<
" not use a marked cut" <<
endl
880 <<
"facei:" << facei <<
endl
881 <<
"face:" <<
f <<
endl
882 <<
"newFace:" << addEdgeCutsToFace(facei)
888 faceCells(cuts, facei, own, nei);
899 faceUptodate[facei] =
true;
907 Pout<<
"meshCutter:" <<
nl
908 <<
" cells split:" << addedCells_.size() <<
nl
909 <<
" faces added:" << addedFaces_.size() <<
nl
910 <<
" points added on edges:" << addedPoints_.size() <<
nl
927 label celli = iter.key();
930 label addedCelli = iter();
934 if (newCelli >= 0 && newAddedCelli >= 0)
939 && (newCelli != celli || newAddedCelli != addedCelli)
942 Pout<<
"meshCutter::topoChange :"
943 <<
" updating addedCell for cell " << celli
944 <<
" from " << addedCelli
945 <<
" to " << newAddedCelli <<
endl;
947 newAddedCells.
insert(newCelli, newAddedCelli);
952 addedCells_.transfer(newAddedCells);
960 label celli = iter.key();
963 label addedFacei = iter();
967 if ((newCelli >= 0) && (newAddedFacei >= 0))
972 && (newCelli != celli || newAddedFacei != addedFacei)
975 Pout<<
"meshCutter::topoChange :"
976 <<
" updating addedFace for cell " << celli
977 <<
" from " << addedFacei
978 <<
" to " << newAddedFacei
981 newAddedFaces.
insert(newCelli, newAddedFacei);
986 addedFaces_.transfer(newAddedFaces);
995 addedPoints_.begin();
996 iter != addedPoints_.end();
1000 const edge&
e = iter.key();
1006 label addedPointi = iter();
1010 if ((newStart >= 0) && (newEnd >= 0) && (newAddedPointi >= 0))
1012 edge newE =
edge(newStart, newEnd);
1017 && (
e != newE || newAddedPointi != addedPointi)
1020 Pout<<
"meshCutter::topoChange :"
1021 <<
" updating addedPoints for edge " <<
e
1022 <<
" from " << addedPointi
1023 <<
" to " << newAddedPointi
1027 newAddedPoints.
insert(newE, newAddedPointi);
1032 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.
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 addCell(const label masterCellID)
Add cell and return new cell index.
label addPoint(const point &, const label masterPointID, const bool inCell)
Add point and return new point index.
label addFace(const face &f, const label own, const label nei, const label masterFaceID, const bool flipFaceFlux, const label patchID)
Add face to cells and return new face index.
#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]